Hybrid quantum-classical computer for variational coupled cluster method

ABSTRACT

A hybrid quantum classical (HQC) computer, which includes both a classical computer component and a quantum computer component, solves linear systems. The HQC decomposes the linear system to be solved into subsystems that are small enough to be solved by the quantum computer component, under control of the classical computer component. The classical computer component synthesizes the outputs of the quantum computer component to generate the complete solution to the linear system.

BACKGROUND

Quantum computers promise to solve industry-critical problems which are otherwise unsolvable. Key application areas include chemistry and materials, bioscience and bioinformatics, and finance. Interest in quantum computing has recently surged, in part, due to a wave of advances in the performance of ready-to-use quantum computers.

Electronic structure problems of strongly correlated systems are among the hardest problems in computational quantum chemistry. They are also of prominent importance because solutions to these problems contain important information about the chemical process under consideration. In the strongly correlated regime, due to entanglement of electronic states, mean field approximations such as Hartree-Fock theory typically fail to capture a significant portion of the ground state energy. In these scenarios, correlated methods such as coupled cluster (CC) and configuration interaction (CI) techniques are used to refine the mean field approximation using physically motivated ansatz wavefunctions. However, both CI and CC have limitations. For example, truncated CI is not size extensive, while full CI methods are not computationally efficient. Coupled cluster methods are able to recover some correlation energy, while its disadvantage is that the coupled cluster energy is not a variational upper bound to the true ground state energy.

Two existing variants of the coupled cluster method, namely unitary coupled cluster (uCC) and variational coupled cluster (vCC), have delivered superior performance to CC, but are too computationally expensive to be implemented on classical computers. It is known that uCC can be implemented efficiently on a quantum computer, thereby making it a viable tool for future quantum chemists once scalable quantum devices have been achieved.

The other variant, however, vCC, has gone largely unexplored. Numerical evidence gleaned from the few instances of vCC implementations show that in the strongly correlated regime at the same level of truncation, vCC could deliver even higher ground state energy than uCC.

What is needed, therefore, are improved techniques for implementing vCC on near-term quantum computers. Such improvements would have a wide variety of applications in science and engineering.

SUMMARY

A hybrid quantum classical (HQC) computer, which includes both a classical computer component and a quantum computer component, implements improved variational coupled cluster (vCC) methods.

In one aspect, a method is performed by a classical computer for implementing, on a quantum computer, a non-unitary operation of the form I+αU, where I is the identity operator, a is a scalar, and U is a unitary operator; the quantum computer having a plurality of qubits, including an ancilla qubit; the classical computer including a processor, a non-transitory computer-readable medium, and computer program instructions stored in the non-transitory computer-readable medium, the computer program instructions being executable by the processor to perform the method. The method includes: (A) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a first quantum circuit W which, when executed by the quantum computer, probabilistically realizes the non-unitary operation by the technique of linear combination of unitaries; (B) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a second quantum circuit, the second quantum circuit comprising a sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1) parametrized by an integer d, wherein S₀=WRW^(†)RW and wherein

$R = \begin{pmatrix} 1 & 0 \\ 0 & e^{i{\pi/3}} \end{pmatrix}$

is a single-qubit rotation applied to the ancilla qubit.

The method may further include: (C) on the quantum computer, executing the first quantum circuit to probabilistically realize the non-unitary operation by the technique of linear combination of unitaries. The method may further include: (D) on the quantum computer, executing the second quantum circuit, comprising executing the sequence of quantum gates S_(d)=S_(d−1)RD_(d−1) ^(†)RS_(d−1). (D) may include, on the quantum computer: (D)(1) applying a single-qubit rotation R_(α) to transform the ancilla qubit into the state

$\left. {\left. {\frac{1}{\sqrt{1 + \alpha}}❘0} \right\rangle + \sqrt{\frac{\alpha}{1 + \alpha}}} \middle| 1 \right\rangle;$

(D)(2) applying controlled-U operator on a subset of the plurality of qubits, not including the ancilla qubit, conditioned on the ancilla qubit being in the state |1

; and (D)(3) applying R_(α) ^(†), on the ancilla qubit.

The method may further include generating and storing, in the non-transitory computer-readable medium, additional computer-readable data that, when executed on the quantum computer, causes the second quantum circuit to execute repeatedly, on the quantum computer, to perform a sequence of operations

${\Pi_{i = 1}^{r}{\Pi_{j = 1}^{k}\left( {1 + {\frac{t\kappa_{j}}{r}P_{j}}} \right)}},$

that approximates the operator e^(tT({right arrow over (κ)})) with t being a scalar and T({right arrow over (κ)})=Σ_(j=1) ^(k)κ_(j)P_(j) being the cluster operator which is a linear combination of operators P_(j) parametrized by weights {right arrow over (κ)}.

The method may further include: (C) performing mean-field approximation to generate and store, in the non-transitory computer-readable medium, computer-readable data representing a description of a third quantum circuit to prepare a reference state; (D) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a parametrized quantum circuit for approximating e^(tT({right arrow over (κ)})); (E) on the quantum computer, executing the third quantum circuit to prepare the reference state; (F) on the quantum computer, applying the parametrized quantum circuit to the reference state to generate the ansatz state |ψ({right arrow over (κ)})

; (G) using the classical computer in cooperation with the quantum computer to measure an energy of the ansatz |ψ({right arrow over (κ)})

; and (H) on the classical computer, iteratively tuning the parameters {right arrow over (κ)} to minimize the energy of the ansatz |ψ({right arrow over (κ)})

.

In another aspect, a system includes: a classical computer the classical computer including a processor, a non-transitory computer-readable medium, and computer program instructions stored in the non-transitory computer-readable medium; a quantum computer comprising a plurality of qubits, including an ancilla qubit; wherein the computer program instructions, when executed by the processor, perform a method for implementing, on the quantum computer, a non-unitary operation of the form I+αU, where I is the identity operator, α is a scalar, and U is a unitary operator. The method includes: (A) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a first quantum circuit W which, when executed by the quantum computer, probabilistically realizes the non-unitary operation by the technique of linear combination of unitaries; (B) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a second quantum circuit, the second quantum circuit comprising a sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1) parametrized by an integer d, wherein S₀=WRW^(†)RW and wherein

$R = \begin{pmatrix} 1 & 0 \\ 0 & e^{i{\pi/3}} \end{pmatrix}$

is a single-qubit rotation applied to the ancilla qubit.

The method may further include: (C) on the quantum computer, executing the first quantum circuit to probabilistically realize the non-unitary operation by the technique of linear combination of unitaries. The method may further include: (D) on the quantum computer, executing the second quantum circuit, comprising executing the sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1). (D) may include, on the quantum computer: (D)(1) applying a single-qubit rotation R_(α) to transform the ancilla qubit into the state

${{\frac{1}{\sqrt{1 + \alpha}}\left. 0 \right\rangle} + {\sqrt{\frac{\alpha}{1 + \alpha}}\left. 1 \right\rangle}};$

(D)(2) applying controlled-U operator on a subset of the plurality of qubits, not including the ancilla qubit, conditioned on the ancilla qubit being in the state |1

; and (D)(3) applying R_(α) ^(\) on the ancilla qubit.

The method may further include generating and storing, in the non-transitory computer-readable medium, additional computer-readable data that, when executed on the quantum computer, causes the second quantum circuit to execute repeatedly, on the quantum computer, to perform a sequence of operations

${\Pi_{i = 1}^{r}{\Pi_{j = 1}^{k}\left( {1 + {\frac{t\kappa_{j}}{r}P_{j}}} \right)}},$

that approximates the operator e^(tT({right arrow over (κ)})) with t being a scalar and T({right arrow over (κ)})=Σ_(j=1) ^(k)κ_(j)P_(j) being the cluster operator which is a linear combination of operators P_(j) parametrized by weights {right arrow over (κ)}.

The method may further include: (C) performing mean-field approximation to generate and store, in the non-transitory computer-readable medium, computer-readable data representing a description of a third quantum circuit to prepare a reference state; (D) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a parametrized quantum circuit for approximating e^(tT({right arrow over (κ)})); (E) on the quantum computer, executing the third quantum circuit to prepare the reference state; (F) on the quantum computer, applying the parametrized quantum circuit to the reference state to generate the ansatz state |ψ({right arrow over (κ)})

; (G) using the classical computer in cooperation with the quantum computer to measure an energy of the ansatz |ψ({right arrow over (κ)})

; and (H) on the classical computer, iteratively tuning the parameters {right arrow over (κ)} to minimize the energy of the ansatz |ψ({right arrow over (κ)})

.

Other features and advantages of various aspects and embodiments of the present invention will become apparent from the following description and from the claims.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 shows a diagram of a system implemented according to one embodiment of the present invention.

FIG. 2A shows a flow chart of a method performed by the system of FIG. 1 according to one embodiment of the present invention.

FIG. 2B shows a diagram illustrating operations typically performed by a computer system which implements quantum annealing.

FIG. 3 shows a diagram of a HQC computer system implemented according to one embodiment of the present invention.

FIG. 4 is a diagram of a hybrid quantum-classical (HQC) computer implemented according to one embodiment of the present invention for realizing a variational quantum eigensolver (VQE) with vCC ansatz.

FIG. 5 is a flowchart of a method performed by the HQC computer of FIG. 4 according to one embodiment of the present invention.

DETAILED DESCRIPTION

Embodiments of the present invention are directed to a hybrid quantum classical (HQC) computer, which includes both a classical computer component and a quantum computer component, and which implements improved variational coupled cluster (vCC) methods. One benefit of the improved vCC methods of embodiments of the present invention is that, unlike uCC methods, they broaden the repertoire of ansatzes within the reach of future quantum chemists with access to scalable fault-tolerant quantum computers.

The potential of quantum computers for simulating quantum systems was first noticed in the 1980s, and has been put on more concrete footing since the late 90s. The ability to directly prepare and measure wavefunctions of molecular systems rather than painstakingly encoding such quantum information into classical computers promises to provide a natural means of performing many quantum chemistry calculations, bypassing the exponential computational requirements of classical algorithms for simulating quantum systems.

Although quantum chemistry is widely understood to be a canonical application for quantum computers, several quantum chemistry problems of interest are known to be hard. A typical example is the QMA-completeness of estimating the ground state energy of a k-local Hamiltonian for k≥2 when the spectral gap scales inverse polynomially in the system size. In spite of such complexity associated with these problems, it is widely believed that physical systems have Hamiltonians that are classically intractable without being hard cases of the QMA class. Furthermore, there are still reasonable classical approximation algorithms that have been developed to study the chemistry of quantum systems.

One attractive class of methods that attack the ground state energy problem chooses a parametrized trial solution to the ground state problem for the given Hamiltonian, and optimizes the parameters by treating the energy expectation value as the objective function to minimize. Ideally, one would hope for this optimization to be performed in a variational framework, which would ensure (by the variational theorem) that the estimate outputted by the algorithm is an upper bound on the true ground state energy of the system. A popular family of trial functions that has been found to give accurate results uses an exponential parametrization and gives rise to the Coupled Cluster (CC) methods. However, performing variational optimization with CC trial functions is computationally intractable on classical computers, generally requiring time exponential in the number of electrons and orbitals, rendering it impractical for most quantum chemistry applications. Several non-variational algorithms have nevertheless been developed and found to be practically useful, despite being expensive to implement on classical computers. These variants require one to solve systems of coupled non-linear equations to obtain an approximation to the optimal parameter values.

These CC methods have played an instrumental role in addressing problems related to quantum systems in the strongly correlated regime. Molecular systems with strong correlation effects typically have wavefunctions of a “multiconfigurational” nature, which means that they have collective excitations coming from several reference wavefunctions chosen for an ab-initio approximation scheme. Both single-reference CC methods that take into account high rank excitations, and multireference CC methods have been developed to tackle such cases. One of the strengths of CC methods is their scaling properties, of especial importance for treating large molecules and solid-state systems. CC methods have been rigorously shown (using (diagrammatic) many body perturbation theory) to be ‘size extensive’, which means that quantities of interest, such as the energy, have the correct scaling in the size, or number of electrons, in the system.

But a significant drawback of traditional CC methods is that the energy expectation obtained by these approximation algorithms is not a variational upper bound on the true energy eigenvalues of the Hamiltonian. This stems from two reasons: the use of a similarity transformed Hamiltonian and an asymmetric expectation value, and the lack of a normalization constraint on the trial wavefunction. Both of these compromises are introduced in order to obtain tractable computational tools. Circumventing these issues to obtain polynomial time classical algorithms is widely believed to be impossible.

Quantum computers, however, can perform computation in some ways that are believed to be beyond classical computing. While classical approaches for electronic structure calculations commonly struggle to maintain and handle the full many-body wavefunctions due to their exponential size, such tasks can be efficiently performed on quantum computers in many cases. A salient example is the recent scheme in variational quantum eigensolvers (VQE) using the unitary coupled cluster (uCC) ansatz, where the parametrized variational trial wavefunction is adjusted by alternating between quantum measurement and classical optimization of parameters such that the energy expectation is minimized. However, another CC variant, the variational coupled cluster method (vCC), has not received much attention in the quantum computing literature due to the non-unitary nature of the operator used to generate the ansatz wavefunction, and intrinsic difficulties in implementing such non-unitary operations on a unitary quantum computer.

Although vCC is computationally expensive on classical computers, there is numerical evidence for small systems to suggest that it yields superior empirical performance over traditional projected CC and uCC. It has been shown that by choosing the model of a Lipkin Hamiltonian, uCC and vCC calculations give very different results for strongly correlated systems. This is also illustrated by benchmarking calculations performed for problems such as computing the potential energy surface of Hydrogen Fluoride, the Nitrogen molecule, and symmetric stretching of the water molecule, where vCC is seen to perform better than uCC, especially at larger bond lengths.

Some VQE methods based on uCC, as well as other variants of CC, have been studied by the quantum computing community in the last few years. In contrast, hybrid quantum-classical algorithms for vCC have not yet been explored, one reason being that the ansatz states are obtained from an initial reference state by applying a non-unitary operator. Embodiments of the present invention include a scheme for efficiently preparing vCC ansatz states on a quantum computer. By integrating this ability with the VQE framework, embodiments of the present invention broaden the set of tools available for treating strongly correlated systems on a quantum computer. Techniques disclosed herein for state preparation combine the two algorithmic ingredients of implementing linear combination of unitaries, and fixed-point oblivious amplitude amplification. Embodiments of the present invention take advantage of the special structure of the coupled cluster operator, writing it as a linear combination of tensor products of single-qubit Pauli matrices in order to provide a quantum algorithm that prepares the CC ansatz wavefunction with asymptotic complexity

. Unlike classical methods for vCC, which incur cost that is exponential in m and n, simply due to the curse of dimensions, the cost of our method is only poly(m,n) for fixed level of excitations

. The energy expectation of these wavefunctions can be estimated to precision ϵ using a number of samples that scales as

${\overset{\sim}{\mathcal{O}}\left( {{{poly}\left( {m,n} \right)} \cdot \frac{1}{\epsilon^{2}}} \right)},$

using standard procedures such as Hamiltonian averaging.

In particular, for vCCSD, the vCC method truncated at double excitations, embodiments of the present invention obtain quantum circuits that prepare the ansatz wavefunctions in time that scales as

(m²n²) or

(m⁴) for a constant orbital filling fraction, matching the asymptotic complexity demonstrated for the uCC method.

The following is organized as follows. First, a theory behind coupled cluster methods of embodiments of the present invention is described. Then an embodiment of a quantum algorithm for preparing parametrized CC ansatz states is described. Then a description is provided of how the ansatz may be used with the VQE for ground state calculations. Then, advantages and drawbacks of embodiments of the present invention are described.

Broadly speaking, the goal of quantum chemistry techniques such as vCC is to find approximations to the ground state energy and wavefunction of a specified Hamiltonian. One example that will be discussed here is the case of a molecular Hamiltonian, although the techniques described herein may be used to treat other kinds of systems as well (e.g., atomic nuclei). Embodiments of the present invention may be applied to the ground state properties of molecular electronic Hamiltonians, simplified by neglecting nuclear motion, which have the general form

$\begin{matrix} {{\hat{H} = {{\sum\limits_{i}{\frac{\nabla_{i}^{2}}{2}{- {\sum\limits_{i}\frac{Z}{\overset{\rightarrow}{r_{l}}}}}}} + {\sum\limits_{ij}\frac{1}{{\overset{\rightarrow}{r_{l}} - \overset{\rightarrow}{r_{J}}}}}}},} & (1) \end{matrix}$

in atomic units, where Z is the nuclear charge, and r_(i) denotes the position of the i^(th) electron.

A quantum system of n electrons may be treated with a basis of m spin orbitals, thereby restricting attention to a finite dimensional subspace of the full domain of the unbounded operator in Eqn. 1. A canonical way of writing down the Hamiltonian restricted to the corresponding finite-dimensional subspace of the Hilbert space of the system is in second quantization, in terms of creation and annihilation operators of the chosen basis

$\begin{matrix} {{\hat{H} = {{\sum\limits_{ij}{h_{ij}{\hat{a}}_{i}^{\dagger}{\hat{a}}_{j}}} + {\frac{1}{2}{\sum\limits_{ijkl}{h_{ijkl}{\hat{a}}_{i}^{\dagger}{\hat{a}}_{j}^{\dagger}{\hat{a}}_{k}{\hat{a}}_{l}}}}}},} & (2) \end{matrix}$

and indeed, any Hamiltonian may be written in this Fock space or occupation number representation. Here the indices range over the spin orbitals and the coefficients h_(ij) and h_(ijkl) are classically precomputed electronic overlap integrals representing the matrix elements of the one particle and two particles terms of the Hamiltonian in the chosen basis. In particular, when the Born-Oppenheimer approximation of neglecting nuclear motion is applied, these coefficients are computed under a particular fixed geometry and may be efficiently performed on a classical computer. Recall that the creation and annihilation operators satisfy fermionic anti-commutation relations

{a _(i) ,a _(j) }:=a _(i) a _(j) +a _(j) a _(i)=0,

{a _(i) ^(†) ,a _(j) ^(†)}=0, {a _(i) ,a _(j) ^(†)}=δ_(ij) I.  (3)

The second quantization approach differs from first quantization in incorporating the anti-symmetry constraint on fermionic states into the definition of the operators rather than the wavefunctions themselves. Furthermore, it has the advantage that the wavefunction may be more efficiently represented in the molecular orbital basis as compared to a cartesian basis or grid chosen for a first quantized approach. As a result, methods using first quantization lead to shallow but wide circuits, whereas those using second quantization lead to deeper circuits but require significantly fewer resources for smaller scale simulations.

To find the ground state of the Hamiltonian in (2), coupled cluster theory introduces an (unnormalized) ansatz wavefunction of the form |ψ

=e^({circumflex over (T)})|Ψ₀

, where |Ψ₀

is a reference state that is suitably chosen for the problem at hand (e.g., as a Slater determinant of Hartree-Fock orbitals). {circumflex over (T)} is called the cluster operator and has the form {circumflex over (T)}={circumflex over (T)}⁽¹⁾+{circumflex over (T)}⁽²⁾+ . . . +{circumflex over (T)}^((n))), where each

is an operator representing the

-th order excitation of electrons from occupied spin orbitals to unoccupied (or virtual) spin orbitals:

$\begin{matrix} {{\hat{T}}^{(\ell)} = {\frac{1}{\left( {l!} \right)^{2}}{\underset{b_{1},\ldots\;,{b_{\ell} \in \;{virtual}}}{\sum\limits_{\alpha_{1},\ldots\;,{\alpha_{\ell} \in \;{occupied}}}}{t_{\alpha_{1}\;\ldots\;\alpha_{\ell}}^{b_{1}\;\ldots\; b_{\ell}}a_{b_{1}}^{\dagger}a_{b_{2}}^{\dagger}\ldots\; a_{b_{\ell}}^{\dagger}a_{\alpha_{\ell}}a_{\alpha_{\ell - 1}}\ldots\mspace{11mu}{a_{\alpha_{1}}.}}}}} & (4) \end{matrix}$

For an n electron system, all terms higher than {circumflex over (T)}^((n)) vanish because there cannot be more than n excitations. In practice, the highest level of excitations

is typically limited to be small. For example, one of the most popular coupled cluster variants is CCSD, which involves only single and double excitations (i.e.,

≤2). The method including triple excitations, CCSDT, or treating such excitations with perturbation theory, CCSD(T), has found widespread praise as the gold standard of quantum chemistry calculations. This is due in large part to the ability of this method to accurately compute nearly 99.5% of the correlation energy (with respect to Full Configuration Interaction calculations) for small to medium sized atoms and molecules, despite its non-variational nature.

The intuition behind using an exponentially parametrized family of trial functions comes from many body perturbation theory and statistical mechanics. The CC wavefunction accounts for static correlations between electrons by assigning non-zero amplitudes to Hilbert space states that represent virtual excitations arising from electronic interactions. Furthermore, the exponential parametrization captures multiple excitations better than truncated Configuration Interaction calculations of the same order, since the exponential operator accounts for clusters of excitations (e.g., a pair of double excitations contributes to the theoretical amplitude of quadruple excitations).

An initial guess for the coupled cluster amplitudes

may be evaluated using many-body perturbation theory, or the Møller-Plesset method. Using such an initial guess to obtain a family of trial wavefunctions, CC theory develops various algorithms to optimize the parameters and converge to an accurate representation of the ground state wavefunction within the ansatz manifold, by treating the energy expectation value of the trial wavefunction (estimated using the known form of the Hamiltonian) as the objective function to be minimized.

Embodiments of the present invention may divide the tasks involved in the vCC method between a quantum processor and a classical processor. At iteration k, the quantum processor takes as input a set of parameter values {right arrow over (t)}^((k))≡{right arrow over (t)} and prepares the CC trial wavefunction |ψ({right arrow over (t)})

using a parametrized circuit U({right arrow over (t)}) that implements the exponential of the cluster operator e^({circumflex over (T)}({right arrow over (t)})) to an initial guess |Ψ₀

(e.g., a Slater determinant of the molecular orbital basis chosen for the vCC approximation), and subsequently evaluates the energy expectation value

E

_({right arrow over (t)}) of this trial wavefunction by performing efficient partial tomography or Hamiltonian averaging. The classical processor runs an optimization routine that computes a new set of parameter values {right arrow over (t)}z^((k+1)), which it passes to the quantum processor. The classical optimization routine treats the energy expectation value passed to it as an objective function to minimize, and uses the quantum processor purely as an oracle for evaluating the function at specified parameter values. This completes one round of the algorithm. Thus, we have the following variational prescription

$\begin{matrix} {\left. {\psi\left( \overset{\rightarrow}{t} \right)} \right\rangle = {e^{T{(\overset{\rightarrow}{t})}}\left. \Psi_{0} \right\rangle}} & (5) \\ {E_{vCC} = {\min\limits_{\overset{\rightarrow}{t}}\frac{\left\langle {{\psi\left( \overset{\rightarrow}{t} \right)}{H}{\psi\left( \overset{\rightarrow}{t} \right)}} \right\rangle}{\left\langle {\psi\left( \overset{\rightarrow}{t} \right)} \middle| {\psi\left( \overset{\rightarrow}{t} \right)} \right\rangle}}} & (6) \end{matrix}$

Since we are performing a variational optimization over the ansatz wavefunction, and since normalization is taken care of by virtue of the state being prepared and measured on the quantum processor, embodiments of the present invention obtain a variational method that will give an upper bound on the ground state energy of the input Hamiltonian. The error in the wavefunction, and hence the energy, depends on the choice of basis wavefunctions and the accuracy of the state preparation step.

The state preparation step is described below. Before the method for this step is described, a description is provided of how the fermionic model is mapped onto a qubit model, using the Jordan-Wigner transformation.

The starting point of quantum algorithms for quantum chemistry problems is to represent the fermionic model on a quantum computer, i.e., a system of qubits and quantum gates. There are several ways to encode fermionic systems into spin systems, most of which use the formalism of second quantization. The fermionic Hamiltonian is isospectrally mapped to a Hamiltonian on the system of spins or qubits, enabling the estimation of energy expectation values indirectly through spin measurements.

Since molecular Hamiltonians are unbounded operators acting on infinite dimensional Hilbert spaces, the numerical methods used to approximately solve molecular structure problems always restrict themselves to a finite dimensional subspace, choosing a suitable basis of fermionic modes that spans this approximation subspace. Fermion to qubit mappings encode fermionic states in Fock space or occupation number space into states of a system of qubits, and the target Hamiltonian H acting on m fermionic modes into a simulator Hamiltonian {tilde over (H)} composed of quantum gates acting on t=f(m) qubits, such that for any isometry V that maps fermionic states to qubit states, we have VH={tilde over (H)}V. Examples of such encodings include the Jordan-Wigner (JW) transformation, the Bravyi-Kitaev (BK) transformation, and the Verstraete-Cirac mapping. Several newer methods have also been proposed, with advantages in the number of qubits required, the number of qubits on which a transformed fermionic creation or annihilation operator acts non-trivially. For simplicity, the description herein focuses on the JW transformation to illustrate embodiments of the vCC algorithm for molecular Hamiltonians. This transformation may be modified suitably to work with most of the encodings mentioned above.

The JW transform maps the fermionic creation and annihilation operators a_(p) ^(†) and a_(q) onto the spin raising and lowering operators on a chain of qubits. These creation and annihilation (or “ladder”) operators in the case of molecular simulations describe the addition or removal of an electron from one of a chosen set of orthonormal molecular basis functions, while the corresponding spin operators describe raising or lowering the spin by one unit. The molecular electrons are indistinguishable, but using the JW transform we can map these onto a string of distinguishable qubits, where qubit j stores the occupation number of orbital j.

The following expressions are obtained for the creation and annihilation operators

$\begin{matrix} {{\left. a_{j}\mapsto{1^{{\otimes j} - 1} \otimes \sigma^{+} \otimes \left( \sigma^{z} \right)^{{\otimes m} - j}} \right.\left. a_{j}^{\dagger}\mapsto{1^{{\otimes j} - 1} \otimes \sigma^{-} \otimes \left( \sigma^{z} \right)^{{\otimes m} - j}} \right.},{where}} & (7) \\ {{\sigma^{+} = {{\left. 0 \right\rangle\left\langle 1 \right.} = \frac{X + {iY}}{2}}}{{\sigma^{-} = {{\left. 1 \right\rangle\left\langle 0 \right.} = \frac{X - {iY}}{2}}},}} & (8) \end{matrix}$

and |0

and |1

are represent an unoccupied and occupied orbital, respectively.

Applying the JW transform to the cluster operator, the JW transform may be written in the form

$\begin{matrix} {{{\hat{T}}^{(i)} = {{\sum\limits_{i_{1},\ldots\;,i_{n}}{t_{i_{1},\ldots\;,i_{n}}{\sigma_{i_{1}}^{\alpha_{1}} \otimes \ldots \otimes \sigma_{i_{n}}^{\alpha_{n}}}}} = {\sum\limits_{j}{t_{j}^{(i)}{\hat{T}}_{j}^{(i)}}}}},} & (9) \end{matrix}$

where the α superscript may take values in {0,z,+,−}, representing the (Pauli) operators {1,Z,o⁺,σ⁻}. The multi-index j=(i₁, . . . ,i₂) may be chosen to subsume all the indices used in the notation introduced for the cluster operator. The number of terms in this linear combination scales polynomially in m and n with the degree being determined by the excitation order at which we truncate.

Using the BK transform, the mapping produces tensor products of Pauli operators that involve only

(log m) terms. However, the JW transform will be described as an example herein for ease of exposition, based on recent evidence that the JW transform is more robust to noise on the quantum processor than the BK transform. Embodiments of the present invention, however, are not limited to use with the JW transform.

To prepare CC trial wavefunctions, we need to implement the exponential of this cluster operator, namely

f{circumflex over (()}T)=e ^({circumflex over (T)}).  (10)

Although the operators {circumflex over (T)}^((i)) do not commute, the individual Pauli terms do. Furthermore, these terms are nilpotent and square to zero. We can take advantage of this fact to write out a simple expression for the exponential of the cluster operator, as we will see below.

Embodiments of the present invention begin by applying the JW transform to the cluster operator, and expressing it as a sum of tensor products of single-qubit spin operators. Writing {circumflex over (T)}=Σ_(ij)t_(j) ^((i)){circumflex over (T)}_(j) ^((i)) as in Eqn. 9, the following two observations hold true: 1) all Pauli terms {circumflex over (T)}_(j) ^((i)) commute with each other, and 2) each term is nilpotent, with ({circumflex over (T)}_(j) ^((i)))²=0. This follows from the fermionic anticommutation relations, noticing that the commutator between any two terms {circumflex over (T)}_(j) ^((i)) and {circumflex over (T)}_(j′) ^((i′)), we know that 1) if they contain at least one a_(i) or a_(j) ^(†) in common, each term in the commutator becomes zero, and 2) if they share no ladder operators corresponding to the same spin orbitals, the commutator is zero using the anticommutation relations—since there are always an even number of operators in a term, the negative signs that arise from the anticommutation relations will cancel each other out.

Thus, we can write

$\begin{matrix} {{e^{T} = {{\prod\limits_{ij}e^{t_{j}^{(i)}{\hat{T}}_{j}^{(i)}}} = {\prod\limits_{ij}\left( {1 + {t_{j}^{(i)}{{\hat{T}}_{j}}^{(i)}}} \right)}}},} & (11) \end{matrix}$

i.e., Trotterisation is exact and the Taylor series truncates at linear order. The total number of terms in this expansion is M(

) (see Eqn. 38 below), and scales as

, which is poly(m,n) for fixed

. Each term in this product may then be implemented separately using fixed-point oblivious amplitude amplification. Embodiments of the present invention may split each term into several slices to ensure that each slice is close to unitary.

Consider each term in the product above: each {circumflex over (T)}_(j) ^((i)) is an m-fold tensor product of single qubit operators (see Eqn. 40 below). Of these, i are σ₊ and i are σ⁻ when the excitation order is i (the rest being Z or 1), making {circumflex over (T)}_(j) ^((i)) non-unitary. To proceed further,

$\sigma_{\pm} = \frac{X \pm {iY}}{2}$

may be multiplied to split the tensor product into a sum of 2^(2i) terms, each with the same coefficient 2^(−2i)t_(j) ^((i)). This results in

$\begin{matrix} {{{1 + {t_{j}^{(i)}{\hat{T}}_{j}^{(i)}}} = {1 + {\frac{t_{j}^{(i)}}{2^{2i}}{\sum\limits_{k = 1}^{2^{2i}}\Theta_{k}^{(i)}}}}},} & (12) \end{matrix}$

where each Θ_(k) ^((i)) is now a unitary, specified as a tensor product of the Pauli matrices {1,X,Y,Z}. This is now in the form of a Linear Combination of Unitaries (LCU). Applying the LCU method to the right-hand side of Eqn. 12, for each {circumflex over (T)}_(j) ^((i)) we get a unitary U_(ij) acting on 2^(2i)+1+m qubits. This unitary has the following action on input states |0² ^(2i) ⁺¹

|ψ

$\begin{matrix} {{{U_{ij}\left. 0^{{2i} + 1} \right\rangle\left. \psi \right\rangle} = {{\frac{1}{1 + t_{j}^{(i)}}\left. 0^{{2i} + 1} \right\rangle\left( {1 + {t_{j}^{(i)}{{\hat{T}}_{j}}^{(i)}}} \right)\left. \psi \right\rangle} + \left. \bot \right\rangle}},} & (13) \end{matrix}$

where 2i+1 ancillas are needed because there are 2^(2i)+1 terms, and the sum of coefficients is 1+t_(j) ^((i)). Some of the Θ_(k) ^((i)) have a factor of i=e^(iπ/2) coming from the σ_(±). This may be dealt with by adding a single ancillary qubit and performing a Y-rotation through π/2, since R_(Y)(π/2)⊗Θ_(k) ^((i))(|+

⊗|ψ

)=e^(iπ/2)|+

⊗Θ_(k) ^((i))|ψ

. By the same argument, phases may always be absorbed to redefine Θ_(k) ^((i)), and henceforth it is assumed that there will always be real positive amplitudes t_(j) ^((i))>0.

Measuring the ancillary qubits, if they are all in the |0

state, results in the normalized state

$\frac{\left( {1 + {t_{j}^{(i)}{\hat{T}}_{j}^{(i)}}} \right)\left. \psi \right\rangle}{{\left( {1 + {t_{j}^{(i)}{\hat{T}}_{j}^{(i)}}} \right)\left. \psi \right\rangle}}$

with success probability

$\begin{matrix} {P_{ij} = {\left. \left( \frac{{\left( {1 + {t_{j}^{(i)}{\overset{\hat{}}{T}}_{j}^{(i)}}} \right)\left. \psi \right\rangle}}{1 + t_{j}^{(i)}} \right. \right)^{2}.}} & (14) \end{matrix}$

It may be assumed that, in practice, the initial guess state |Ψ₀

has non-zero overlap with the target ground state. Under this assumption, the numerator is always non-zero.

Embodiments of the present invention may successively apply the U_(ij)s corresponding to each of the terms in (11), and relegate post-selection to the last step. However, proceeding in this fashion leads to an exponentially small success probability. To sidestep this issue, embodiments of the present invention may perform each term using fixed-point amplitude amplification (FPOAA).

Each term e^(t) ^(j) ^((i)) ^({circumflex over (T)}) ^(j) ^((i)) =1+t_(j) ^((i)){circumflex over (T)}_(j) ^((i))is realized using the corresponding U_(ij), with success probability bounded below by the worst-case value of 1/(1+τ)², with

$\tau:={\max\limits_{ij}{t_{j}^{(i)}.}}$

Since the probabilities for different terms multiply, the total runtime in the Repeat-Until-Success (RUS) framework will also be multiplicative and can scale exponentially as the number of terms in Eqn. 11, without using amplitude amplification.

Instead, embodiments of the present invention may split each term into a number of segments by the prescription e^(tA)=(e^(ηA))^(t/η), where η∈(0,1) is a small parameter and r=t/η is the number of segments (for notational convenience we represent a generic term t_(j) ^((i)){circumflex over (T)}_(j) ^((i)) here by tA). Embodiments of the present invention may adjust r such that each segment is

(η)-close to a unitary in the spectral norm, satisfying the conditions required for oblivious amplitude amplification, and the parameter η will enter into the error analysis for the wavefunction. Now using the FPOAA technique, embodiments of the present invention may implement each segment with success probability 1−ε using

(log(1/ε)/log(1/η)) queries. So, for t/η segments, embodiments of the present invention may use

$\mathcal{O}\left( {\frac{t}{\eta}\frac{\log\left( {1/ɛ} \right)}{\log\left( {1/\eta} \right)}} \right)$

calls to the LCU unitary for A. The RUS runtime for a single term in Eqn. 11 then scales as

$\begin{matrix} {{\mathcal{O}\left( {\left( {1 - ɛ} \right)^{\frac{t}{\eta}} \cdot \frac{t}{\eta} \cdot \frac{\log\;\left( {1/ɛ} \right)}{\log\;\left( {1/\eta} \right)}} \right)}.} & (15) \end{matrix}$

If ϵ=η/t, then (1−ε)^(−t/η)=O(1) for η/t<<1.

Consider implementing each term successively using the method outlined above. Let

$\tau = {\max\limits_{ij}t_{j}^{(i)}}$

as before. To get the overall complexity, note that the probabilities multiply across all segments and all terms e^(t) ^(j) ^((i)) ^({circumflex over (T)}) ^(j) ^((i)) , since we need to discard and re-initialize the LCU ancillary qubits for each segment of each term. Thus, if the probability for each segment is amplified by FPOAA to 1−ϵ, the overall success probability becomes

p≥(1−

,

where r=τ/η is the maximum number of segments required for any of the individual terms, and M(

) is the total number of such terms (see Eqn. 38 below). Since FPOAA converges monotonically to the target state, the number of segments r may be chosen for the term having the largest coefficient r=t_(j) ^((i)), and fix the same number of segments for all the terms. Now, if

${\epsilon = \frac{1}{{rM}(\ell)}},$

then

p≥(1−1/rM(

)

.

This converges up towards 1/e>0.36 as rM(

)→∞, and already for rM(

)=10, p>0.35. By choosing ϵ=CrM(

) for some positive integer C>1, the success probability can be made

${p \approx \sqrt[C]{1/e}},$

and for C=2, this gives p≥0.6. Thus, embodiments of the present invention yield a constant success probability greater than ½, and the total number of uses of all the different LCU unitaries U_(ij) for the M(

) terms in Eqn. 11 adds up to

$\begin{matrix} {{\mathcal{O}\left( {{{rM}(\ell)} \cdot \frac{\log\;\left( {{rM}(\ell)} \right)}{\log\;\left( {\tau/r} \right)}} \right)}.} & (16) \end{matrix}$

To choose r, or the equivalent η=τ/r, note that the analysis herein is valid if η<0.5. Thus, for a given set of parameter values of {right arrow over (t)}, r=Cτ may be chosen for some constant C>2 and τ=max_(ij) t_(j) ^((i)) as defined previously is the max-norm of the vector of parameters. Henceforth, the term log(τ/r) may be dropped from asymptotic expressions.

The above expression may be used to get a precise gate count for the case of finite m,n,

. For asymptotics, note that

, so that an asymptotic scaling is given by

,  (17)

hiding the logarithmic factor log(n(m−n)) in the

notation, and suppressing the weak log τ dependence.

Finally, the gate complexity of the circuit can be obtained by noting that the {circumflex over (T)}_(j) ^((i)) term of excitation order i consists of 2^(i) non-trivial gates, which are one of the three Pauli operators X, Y, and Z. Thus, a simple upper bound on the gate complexity is given by the expression in Eqn. 16 multiplied by the factor

. The actual circuit may be simpler than the gate count indicates, since Z and X gates are easy to implement, and it may be possible to combine terms and reuse circuit segments.

There are several ways to choose single and multireference initial guess states |Ψ₀

using Hartree-Fock or MCSCF calculations as a starting point. For the particular implementation described in FIG. 4, after canonical basis functions are produced using Hartree-Fock the reference state is

$\left. {\underset{\underset{n}{︸}}{11{\cdots 1}}00{\cdots 0}} \right\rangle$

recalling n being the number of electrons. Embodiments of the present invention may easily prepare this on a quantum computer.

Ground State Calculations Using the VQE.

The first step in using the VQE framework is preparing the ansatz wavefunction on the quantum processor. This may be followed by an efficient estimation of the average energy of this trial wavefunction. Embodiments of the present invention may, for example, use a Hamiltonian averaging procedure to estimate this energy.

The Hamiltonian is mapped by the same isospectral mapping to a qubit operator. For example, this may be in the form of a sum of m-fold tensor products of Pauli operators in the case of the JW transformation

$\begin{matrix} {{H = {\sum\limits_{i}{h_{i}O_{i}}}},} & (18) \end{matrix}$

corresponding to the Hamiltonian in Eqn. 2, where the coefficients h_(i) are the precomputed overlap integrals (matrix elements of the Hamiltonian in the basis chosen for the vCC method). The number of terms in this sum scales as

(m⁴) for molecular Hamiltonians, which contain at most two-body interactions. The average energy of the trial wavefunction |ψ({right arrow over (t)})

is given by

$\begin{matrix} {{{\left\langle H \right\rangle_{\overset{\rightarrow}{t}}:} = {\left\langle {{\psi\left( \overset{\rightarrow}{t} \right)}{H}{\psi\left( \overset{\rightarrow}{t} \right)}} \right\rangle = {\sum\limits_{i}{h_{i}\left\langle O_{i} \right\rangle_{\overset{\rightarrow}{t}}}}}},} & (19) \end{matrix}$

which may be determined purely by making single qubit spin measurements, each of which have only two outcomes ±1. This in turn means that the variance of each of the observables O_(i) is bounded by 1. If we make M_(i) measurements to determine the quantity

O_(i)

, the precision to which we can approximate the value of the corresponding term in the Hamiltonian is given by the empirical variance of the mean we estimate, that is

$\begin{matrix} {\epsilon_{i}^{2} = {\frac{{h}_{i}^{2}\;{{Var}\left( \left\langle O_{i} \right\rangle \right)}}{M_{i}} \leq {\frac{{h}_{i}^{2}}{M_{i}}.}}} & (20) \end{matrix}$

To obtain the total energy to a precision ϵ² using a total of M=Σ_(i)M_(i) measurements (or equivalently, samples), embodiments of the present invention may choose

$\begin{matrix} {{M_{i} = {\frac{h_{i}}{\underset{j}{\Sigma}{h_{j}}}M}},} & (21) \end{matrix}$

and since the total error in the energy estimated using these measurements is given by the sum of the variances, Σ_(i)e_(i) ², this leads to

$\begin{matrix} {{M \approx \frac{\left( {\Sigma_{i}{h_{i}}} \right)^{2}}{\epsilon^{2}}},} & (22) \end{matrix}$

so that the number of measurements, and hence samples of the trial wavefunction, scales quadratically with the inverse precision to which we would like to estimate the energy. While it does require a fairly large number of measurements, this method makes it possible to avoid quantum phase estimation, which requires long coherence times and the use of

(log 1/ϵ) ancillary qubits.

Thus, the net complexity of the part of the embodiments of the present invention that run on the quantum co-processor, including the number of measurements required for Hamiltonian averaging to find the energy expectation value of trial wavefunctions, has an asymptotic scaling of

${\overset{\sim}{\mathcal{O}}\left( {{\frac{1}{ɛ^{2}} \cdot \tau}\;\ell\; m^{\ell\;}n^{\ell\;}} \right)},$

where the factor (Σ_(i)|h_(i)|)² has been hidden, which is related to the norm of the input Hamiltonian, and scales as

(m⁴) in the worst case.

The next step in a round of variational quantum optimization is to perform classical optimization over the parameter space of {right arrow over (t)} using the average energy of the trial wavefunction as the objective to minimize. A variety of methods, ranging from stochastic gradient descent to gradient-free direct search methods such as the Nelder-Mead simplex algorithm, may be used for this purpose. The cost of optimization on the classical processor will depend on the dimension of the parameter space, which in the worst-case scales asymptotically as

. Thus, for vCCSD (vCC truncated to single and double excitations), the cost depends on

(m²n²). In practice, the parameter space may, for example, be pruned by using preprocessing using Møller-Plesset or many body perturbation theory calculations to estimate the cluster amplitudes and remove those terms from the cluster operator that have amplitude below some threshold.

In addition to the number of parameters to optimize, the number of rounds required by the optimization routine to achieve convergence will depend on the quality of the initial guess or reference state. In particular, one can expect it to depend inverse polynomially on the magnitude of the overlap β=

tar|Ψ₀

between the reference state and the true target state, which in the simplest case is the true ground state. This may, for example, be tuned to some extent by choosing suitable reference states, using Hartree-Fock, MP2, MCSCF or other ab-initio techniques. Indeed, unless the reference state has |β|≥1/poly(m,n), the optimization cannot be performed in poly(m,n) time.

As described above, embodiments of the present invention include a quantum algorithm for preparing the exponentially parametrized coupled cluster family of ansatz states from variational quantum chemistry with asymptotic gate complexity

for a constant excitation order

. Using a Hamiltonian averaging procedure, the energy expectation value of the trial functions thus prepared may be estimated to precision ϵ using only single qubit spin measurements, and a total of

(Σ_(i)|h_(i)|²/ϵ²)=

(m⁴/ϵ²) samples. Here the h_(i) are the matrix elements of the fermionic Hamiltonian in the molecular orbital basis chosen for the vCC approximation, and the loose bound of

(m⁴) is obtained simply by counting the number of terms in the second-quantised Hamiltonian. This is exponentially better than known classical methods of estimating the energy expectation value of the trial wavefunction, which suffer from the usual problem of having to deal with vectors of dimension exponential in m, particularly in the strongly correlated regime where the target state is expected to be highly entangled and delocalised. These results match the performance of other known CC methods.

In the special case of vCC truncated to single and double excitations, vCCSD, the entire quantum algorithm disclosed herein for evaluating the energy expectation value of the trial wavefunction may have a complexity that scales as

(m²n²·1/ϵ²).

Embodiments disclosed herein may, for example, be modified in a variety of ways, such as by using multireference initial guess states, and extended coupled cluster techniques that use modified cluster operators and include de-excitation terms.

The Linear Combination of Unitaries (LCU) method is a technique that can be used to probabilistically implement any operator that is decomposed into a linear combination of unitaries, each of which can be implemented.

Given an operator A=Σ_(i=1) ^(N) t_(i) Θ_(i), where the Θ_(i) are unitary, A may be implemented, for example, using the LCU method as follows. First, two unitaries may be constructed, one that performs the controlled operation

$\begin{matrix} {{{{SELECT}\mspace{14mu}(A)} = {\sum\limits_{j = 0}^{N - 1}{\left. j \right\rangle{\left\langle j \right. \otimes \Theta_{j}}}}},} & (23) \end{matrix}$

and one that prepares a control-state, performing the map

$\begin{matrix} {{{V\left. 0^{m} \right\rangle} = {\frac{1}{\sqrt{\tau}}{\sum\limits_{j = 0}^{N - 1}{\sqrt{t_{j}}\left. j \right\rangle}}}},} & (24) \end{matrix}$

where τ=Σ_(j)|t_(j)|, and the number of ancillary qubits required is m=┌log N┐. Then the LCU operator may be defined as

W=(V ^(†)⊗1)SELECT(A)(V⊗1),  (25)

where the identity part in the first and third terms acts on the system register. The action of W on a state |Ψ

=|0^(m)

|ψ

is given by

$\begin{matrix} {{{W\left. 0^{m} \right\rangle\left. \psi \right\rangle} = {{\frac{1}{\tau}\left. 0^{m} \right\rangle A\left. \Psi \right\rangle} + \left. \bot \right\rangle}},} & (26) \end{matrix}$

where (|0^(m)

0^(m)|⊗1)|⊥

=0. If we postselect on the ancillary register being in the state |0^(m)

, the output state will be

$\begin{matrix} {{\left. \psi_{out} \right\rangle = \frac{A\left. \psi \right\rangle}{{A\left. \psi \right\rangle}}},} & (27) \end{matrix}$

and the probability of success in obtaining this state will be

$\begin{matrix} {{p = {{\frac{{A❘\psi}}{\tau}}^{2} \geq {\frac{A^{- 1}}{\tau}}^{2}}},} & (28) \end{matrix}$

assuming A is not singular, where by ∥A⁻¹∥ we mean the least eigenvalue (singular value) of A on the domain of interest. Depending on the setup, different versions of amplitude amplification may be used to boost the success probability using

$\begin{matrix} {{\mathcal{O}\left( \frac{1}{\sqrt{p}} \right)} = {\mathcal{O}\left( {\frac{\tau}{A^{- 1}}} \right)}} & (29) \end{matrix}$

repetitions of the LCU operator W. For an efficient implementation, the value of the sum of coefficients τ must not be too large, and ∥A⁻¹∥ must not be too small.

Next, examples of techniques for preparing the control state for the LCU will be described. To implement T:=T_(j) ^((i)) given by:

$T = {1 + {\frac{t}{2^{2i}}{\sum\limits_{l = 0}^{2^{2i} - 1}\Theta_{j}^{(i)}}}}$

using the LCU method, we need a control-state preparation unitary that performs the map:

${{v\left. 0 \right\rangle} = {\frac{1}{\sqrt{1 + t}}\left( {\left. 0 \right\rangle + {\sum\limits_{l = 1}^{2^{2i} - 1}{\sqrt{\alpha}\left. l \right\rangle}}} \right)}},$

where t:=t_(j) ^((i)) and

$\alpha = {\frac{t}{2^{2i}}.}$

Such a unitary can be constructed using known techniques. Consider the probability distribution p=(p₀, . . . ,p₂ _(2i+1) ) with:

$p_{0} = \frac{1}{1 + t}$ $P_{i} = {{\frac{\alpha}{1 + t}1} \leq i \leq 2^{2i}}$ P_(k) = 02^(2i) + 1 ≤ k ≤ 2^(2i + 1).

This distribution may readily be integrated (in this case, summed up) to obtain the cumulative distribution function, which is required in creating the quantum state corresponding to this p. If 0, . . . ,2^(2i+1) is partitioned into 2^(k) segments of equal size, numbered 0, . . . ,2^(k)−1, and consider the conditional probability concentrated in the left half of the l^(th) segment, this results in:

${f_{k}(0)} = {{\frac{1}{p_{k - 1}^{0}} \cdot \frac{1}{1 + t}}\left( {1 + {\left( {2^{{2i} - k} - 1} \right)\alpha}} \right)}$ ${f_{k}(l)} = {{{\frac{1}{p_{k - 1}^{l}} \cdot \frac{1}{1 + t}}\left( {2^{{2i} - k}\alpha} \right)1} \leq l \leq {2^{k - 1} + 1}}$ f_(k)(l) = 02^(k − 1) + 2 ≤ l ≤ 2^(k),

where the normalization factors are the segment probabilities from the previous step. Then, embodiments of the present invention may perform the following procedure. The following single qubit state may be prepared:

|ψ₁

=(√{square root over (f(0))}|0

+√{square root over (f(1))}|1

),

corresponding to dividing the set 0, . . . ,2^(2i+1) into two equal segments. In this case,

${f(0)} = {{1 - {\frac{\alpha}{1 + t}\mspace{14mu}{and}\mspace{14mu}{f(1)}}} = {\frac{\alpha}{1 + t}.}}$

Embodiments of the present invention may then adjoin an ancilla qubit in the state |0

and, by applying conditional rotations, prepare a state corresponding to 2²=4 segments. The following conditional single-qubit rotation may be applied:

|0

0|⊗R _(y)(2 cos⁻¹ √{square root over (f ₁(0))})+|1

1|⊗R _(y)(2 cos⁻¹ √{square root over (f ₁(1))}),

where

R _(y)(2 cos⁻¹ √{square root over (f _(k)(l))})|0

=√{square root over (f _(k)(l))}|0

+√{square root over (1−f _(k)(l))}|1

.

At step k, when there is the state corresponding to 2^(k) segments:

${\left. \psi_{k} \right\rangle = {\sum\limits_{l = 0}^{2^{k} - 1}{\sqrt{p_{k}^{l}}\left. l \right\rangle}}},$

where p_(k) ^(l) is the probability contained in the l^(th) segment, embodiments of the present invention may apply the following conditional rotation to go to step k+1, i.e., to further subdivide into 2^(k+1) segments:

$\sum\limits_{l = 0}^{2^{i} - 1}{\left. l \right\rangle{\left\langle l \right. \otimes {{R_{y}\left( {2\cos\limits^{- 1}\sqrt{f_{k}(l)}} \right)}.}}}$

Thus, for excitation order i, exactly 2i+1 ancillary qubits may be used, and 2^(2i)+1 non-trivial controlled single-qubit rotations. If embodiments of the present invention truncate at excitation order

, then the total number of ancillas and controlled rotations is

and

respectively.

For CCSD and CCSDT, this corresponds to 5 or 7 ancillae, and 16 or 64 single-qubit rotations per state preparation.

Fixed-point oblivious amplitude amplification (FPOAA), which has recently been investigated, is a method that combines the features of oblivious amplitude amplification and fixed-point quantum search. The advantage of the former is that, unlike traditional amplitude amplification, one does not need to be able to reflect about the entire target state. The fixed-point search algorithm, on the other hand, has the advantage of monotonously converging to the target state, avoiding the overshooting problem of traditional amplitude amplification.

Suppose there is a unitary W which implements an operation

W|0^(m)

|ψ

=α|0^(m)

V|ψ

+β|Φ ^(⊥)

where V is a linear combination of unitaries that is also a unitary. Such a W preserves a 2-dimensional space spanned by |Φ

=|0^(m)

V|ψ

and |Φ^(⊥)

. Let Π=|0^(m)

0^(m)|0^(m)⊗1 be a projection onto the ancillary part of the target state, and consider the operator R=I−(1−e^(iπ/3))Π that applies a selective phase shift. Let S=WRW^(†)RW. Then if |α|²=1−ϵ, one can show that music: ∥ΠS|0^(m)

|ψ∥²=1−ϵ³. This probability amplification method can also be constructed recursively, as follows

S _(d) =S _(d−1) RS _(d−1) ^(†) RS _(d−1) , S ₀ =S,  (30)

where recursion to depth d decreases the failure probability from ϵ to ϵ³ ^(d) , requiring

$\frac{3^{d} - 1}{2}$

calls to W.

Hence to suppress the failure probability to at most E from an initial value of ε₀, one needs O(log ε⁻¹/log ε₀ ⁻¹) calls to W.

Embodiments of the present invention may implement a non-unitary operation {tilde over (V)}=I+tT using the operation {tilde over (W)}=U_(ij) in (13). For small t, {tilde over (V)} is close to unitary and the same procedure as above works, with a suitably modified error analysis. Hence

$\begin{matrix} {{\overset{\sim}{W}\left. 0^{m} \right\rangle\left. \psi \right\rangle} = {{\frac{1}{1 + t}\left. 0^{m} \right\rangle\overset{\sim}{V}\left. \psi \right\rangle} + \left. {\overset{\sim}{\Phi}}^{\bot} \right\rangle}} & (31) \end{matrix}$

where |{tilde over (Φ)}^(⊥)

is some unnormalized state. Let {tilde over (S)}={tilde over (W)}R{tilde over (W)}^(†)R{tilde over (W)}. Noticing that the success probability is

$p = {\frac{{\overset{\sim}{V}\left. \psi \right\rangle}}{1 + t}}^{2}$

and defining

${\left. \Phi \right\rangle = {{\frac{\left. 0^{m} \right\rangle\overset{\sim}{V}\left. \psi \right\rangle}{{\left. 0^{m} \right\rangle\overset{\sim}{V}\left. \psi \right\rangle}}\mspace{14mu}{and}\mspace{14mu}\left. \Phi^{\bot} \right\rangle} = \frac{\left. {\overset{\sim}{\Phi}}^{\bot} \right\rangle}{\left. {\overset{\sim}{\Phi}}^{\bot} \right\rangle }}},$

embodiments of the present invention may take advantage of unitarity to rewrite Eqn. 31 with the states on the right-hand side normalized as

{tilde over (W)}|0^(m)

|ψ

=√{square root over (p)}|Φ

+√{square root over (1−p)}|Φ^(⊥)

  (32)

The result of this, in terms of the normalized states, is

$\begin{matrix} {{\overset{\sim}{S}\left. 0^{m} \right\rangle\left. \psi \right\rangle} = {\overset{\sim}{W}R{\overset{\sim}{W}}^{\dagger}{R\left( {{\sqrt{p}\left. \Phi \right\rangle} + {\sqrt{1 - p}\left. \Phi^{\bot} \right\rangle}} \right)}}} \\ {= {\overset{\sim}{W}R{{\overset{\sim}{W}}^{\dagger}\left( {{e^{i{\pi/3}}\sqrt{p}\left. \Phi \right\rangle} + {\sqrt{1 - p}\left. \Phi^{\bot} \right\rangle}} \right)}}} \\ {= {\overset{\sim}{W}{R\left\lbrack {{\left( {e^{i{\pi/3}} - 1} \right)\sqrt{p}{\overset{\sim}{W}}^{\dagger}\left. \Phi \right\rangle} + {\left. 0^{m} \right\rangle\left. \psi \right\rangle}} \right\rbrack}}} \\ {= {{\overset{\sim}{W}\left( {I - {\left( {1 - e^{i{\pi/3}}} \right)\prod}} \right)}\left\lbrack {{\left( {e^{i{\pi/3}} - 1} \right)\sqrt{p}{\overset{\sim}{W}}^{\dagger}\left. \Phi \right\rangle} + {\left. 0^{m} \right\rangle\left. \psi \right\rangle}} \right\rbrack}} \\ {= {{\left( {e^{i{\pi/3}} - 1} \right)\sqrt{p}\left. \Phi \right\rangle} + {e^{i{\pi/3}}\overset{\sim}{W}\left. 0^{m} \right\rangle\left. \psi \right\rangle} + {\sqrt{p}\left( {e^{i{\pi/3}} - 1} \right)^{2}\overset{\sim}{W}{\prod{{\overset{\sim}{W}}^{\dagger}\left. \Phi \right\rangle}}}}} \end{matrix}$

The third step above uses Eqn. 32 to write |0^(m)

|ψ−√{square root over (p)}{tilde over (W)}^(†)|Φ

=√{square root over (1−p)}{tilde over (W)}^(†)|Φ^(⊥)

. The last term on the last line may be simplified to

{tilde over (W)}Π{tilde over (W)} ^(†)|0^(m)

{tilde over (V)}|ψ

=p|Φ

+√{square root over (p(1−p))}|Φ^(⊥)

.

This results in

{tilde over (S)}|0^(m)

|ψ

=|Φ

(√{square root over (p)}(2e ^(iπ/3)−1)+(√{square root over (p)})³(e ^(iπ/3)−1)²)+|Φ^(⊥)

(√{square root over (1−p)}e ^(iπ/3) +p√{square root over (1−p)}(e ^(iπ/3)−1)²).  (34)

The failure probability, or probability of obtaining a measurement result in the non-target subspace, is given by

p′=|√{square root over (1−p)}e ^(iπ/3) +p√{square root over (1−p)}(e ^(iπ/3)−1)²|²=(1−p)³,  (35)

which shows that the failure probability falls from 1−p to (1−p)³.

Furthermore,

$p \geq \frac{1}{\left( {1 + t} \right)^{2}} \geq {1 - {2t}}$

for t<½. Therefore, the initial failure probability ε₀≤2t, and the total number of calls to W required to suppress the failure probability down to E scales as O(log(1/ε)/log(1/t)).

For an n-electron system, and a basis

=

_(o)∪

_(u) consisting of a total of m single particle orbitals, of which |

_(o)|=n are occupied and |

_(u)|=m−n are unoccupied, the cluster operator term of excitation order j has the form

$\begin{matrix} {{\hat{T}}^{(j)} = {\frac{1}{\left( {j!} \right)^{2}}{\sum\limits_{i_{1},\ldots\mspace{14mu},{i_{j} \in \mathcal{B}_{o}}}{\sum\limits_{\alpha_{1},\ldots\mspace{14mu},{\alpha_{j} \in \mathcal{B}_{u}}}{t_{\alpha_{1}\alpha_{2}\mspace{14mu}\ldots\mspace{14mu}\alpha_{j}}^{i_{1}i_{2}\mspace{14mu}\ldots\mspace{14mu} i_{j}}a_{\alpha_{1}}^{\dagger}a_{\alpha_{2}}^{\dagger}\mspace{14mu}\ldots\mspace{14mu} a_{\alpha_{j}}^{\dagger}a_{i_{j}}a_{i_{j - 1}}\mspace{14mu}\ldots\mspace{14mu}{a_{i_{1}}.}}}}}} & (36) \end{matrix}$

Since there are j creation and annihilation operators each, corresponding to the occupied and unoccupied single particle basis states, the total number of terms in this sum is given by (neglecting constant factors that may arise due to normal ordering)

$\begin{matrix} {M_{j}:={{{Num}\left( {\hat{T}}^{(j)} \right)} = {{\begin{pmatrix} n \\ j \end{pmatrix}\begin{pmatrix} {m - n} \\ j \end{pmatrix}} = {{\mathcal{O}\left( {n^{j}\left( {m - n} \right)}^{j} \right)}.}}}} & (37) \end{matrix}$

If these are summed up for excitations up to order

, the result is

$\begin{matrix} {{{M(\ell)}:={{\sum\limits_{i = 0}^{\ell}{\begin{pmatrix} n \\ i \end{pmatrix}\begin{pmatrix} {m - n} \\ i \end{pmatrix}}} = {\frac{\Gamma\left( {m + 1} \right)}{{\Gamma\left( {n + 1} \right)}{\Gamma\left( {m - n + 1} \right)}} - {\begin{pmatrix} n \\ {\ell + 1} \end{pmatrix}\begin{pmatrix} {m - n} \\ {\ell + 1} \end{pmatrix}\mspace{34mu}{{{}_{}^{}{}_{}^{}}\left( {1,{\ell - n + 1},{{\ell - m + n + 1};{\ell + 2}},{{\ell + 2};1}} \right)}}}}},} & (38) \end{matrix}$

where ₃F₂ is a generalised hypergeometric function.

Although asymptotically the number of terms scales as

, the polynomial that is hiding under the highest degree term has a large number of terms and potentially large coefficients. For practical purposes, it would also be useful to know an exact gate count.

Furthermore, note that the number of terms also equals the number of independent parameters that have to be optimized over by the classical optimization routine in the VQE procedure.

One ingredient that embodiments of the present invention may use in preparing CC ansatz wavefunctions using the LCU method is the expansion of the JW transformed cluster operator as a linear combination of tensor products of single qubit Pauli operators, or “Pauli strings.”

Consider {circumflex over (T)}=Σ_(i){circumflex over (T)}^((i))=Σ_(ij){circumflex over (T)}_(j) ^((i)), in which the upper index represents excitation order. The lower index runs from 1, . . . ,M_(i) for M_(i) from Eqn. 37 and labels an ordered pair of multi-indices ({right arrow over (b)},{right arrow over (α)}) which are both i-tuplets ordered so that b₁<b₂< . . . <b_(i)<α₁<α₂< . . . <α_(i), where 1≤b_(k)≤n and n<α_(k)≤m for all 1≤k≤i. The operator {circumflex over (T)}^((i)) may be written as

$\begin{matrix} {{\hat{T}}^{(i)} = {{\sum\limits_{j = 1}^{M_{i}}T_{j}^{(i)}} = {\sum\limits_{{b_{1} < b_{2} < \ldots < b_{i}}{\forall{{kb}_{k} \in \mathcal{B}_{o}}}}{\sum\limits_{{\alpha_{1} < \alpha_{2} < \ldots < \alpha_{i}}{\forall{{k\alpha}_{k} \in \mathcal{B}_{u}}}}{t_{\overset{\rightarrow}{\alpha}}^{\overset{\rightarrow}{b}}a_{b_{1}}^{\dagger}a_{b_{2}}^{\dagger}\mspace{14mu}\ldots\mspace{14mu} a_{b_{i}}^{\dagger}a_{\alpha_{i}}a_{\alpha_{i - 1}}\mspace{14mu}\ldots\mspace{14mu}{a_{\alpha_{1}}.}}}}}} & (39) \end{matrix}$

Now applying the JW transform (see Eqn. 7) for each of the a,a^(†) terms and using Zσ_(±)=±σ_(±) etc., the result is

$\begin{matrix} {{\hat{T}}_{j}^{(i)} = {1^{{\otimes b_{1}} - 1} \otimes \sigma_{-} \otimes Z^{{\otimes b_{2}} - b_{1} - 1} \otimes \left( {- \sigma_{-}} \right) \otimes 1^{{\otimes b_{3}} - b_{2} - 1} \otimes \mspace{14mu}\ldots \otimes \left( {Z^{i - 1}\sigma_{-}} \right) \otimes \left\{ \begin{matrix} {{{Z^{{\otimes \alpha_{1}} - b_{i} - 1} \otimes \sigma_{+} \otimes 1^{{\otimes \alpha_{2}} - \alpha_{1} - 1} \otimes \sigma_{+} \otimes \mspace{14mu}\ldots\mspace{14mu} \otimes \sigma_{+}}1^{m - \alpha_{i}}},} & i \\ {{{1^{{\otimes \alpha_{1}} - b_{i} - 1} \otimes \sigma_{+} \otimes Z^{{\otimes \alpha_{2}} - \alpha_{1} - 1} \otimes \sigma_{+} \otimes \mspace{14mu}\ldots\mspace{14mu} \otimes \sigma_{+}}1^{m - \alpha_{i}}},} & i \end{matrix} \right.}} & (40) \end{matrix}$

Note that the Pauli string {circumflex over (T)}_(j) ^((i)) may further be indexed by an m-bit binary string—simply by representing both α_(±) by 1 and both 1 and Z by 0. Then it is known that:

if the string has Hamming weight 2i, then i of the 1 s must appear within the first n places, and the remaining i in the last m−n places;

the first i instances of 1 correspond to σ⁻ (i.e., to the creation operators a_(b) _(k) ^(†)) and the subsequent i 1 s correspond to σ₊ (i.e., to the annihilation operators a_(a) _(k) );

there can be at most 2n positions with a 1.

Furthermore, the mapping f:⊗_(i=1) ^(m) σ_(i)→∪_(i=1) ^(n) S_(2i)⊂{0,1}^(m) of strings of the form of Eqn. 40 into 7n-bit binary strings of even weight up to 2n gives a way to order the terms within each excitation order {circumflex over (T)}^((i)), or to order all terms within {circumflex over (T)} as a whole, whichever we choose to do. That is, a natural family of indexing functions is obtained by setting

:f({circumflex over (T)}_(j) ^((i)))=f₁ . . . f_(m)

k∈[M_(i)] if the nt-bit string f of weight 2i represents in binary the k^(th) smallest integer in the list of weight 2i strings allowed under the rules described above.

Alternatively, a function

may be obtained that simply sorts all the strings allowed by the above rules (across all weights) in ascending order, assigning a string the integer corresponding to its position in the sorted list. These indexing functions and their inverses are efficiently computable in time

(poly(m,n)). For example, by using the [n][m−n] structure of the Hamming weight 2i strings, the operator can be reconstructed from the corresponding bit string on the fly.

A HQC 400 according to one embodiment of the present invention is shown in FIG. 4. In general, the HQC 400 realizes a variational quantum eigensolver (VQE) with vCC ansatz, as will now be described in more detail. The HQC 400 includes both a classical computer 402 and a quantum computer 104.

For a quantum mechanical problem involving n interacting electrons, one typically considers an electronic structure Hamiltonian H of the form:

$\begin{matrix} {H = {{\sum\limits_{i = 1}^{n}T_{i}} + U_{i} + {\sum\limits_{i,j,{i \neq j}}V_{ij}}}} & (41) \end{matrix}$

where T_(i)=½∇_(i) and U_(i) are the kinetic and potential energies of each electron, respectively, and

$V_{ij} = \frac{1}{{r_{i} - r_{j}}}$

is the Coulomb interaction energy between electron i at coordinate r_(i) and electron j at coordinate r_(j). Hence, the electronic structure Hamiltonian H is a continuous operator acting on the coordinate space r. The spectrum of H in most cases is beyond analytical solution. One general method for approximating the spectrum of H is to introduce a set of orthonormal basis functions

$\begin{matrix} {\left\{ {{\psi_{i,s}(r)},{i = 1},\ldots\mspace{14mu},\frac{N}{2},{s = {- 1}},{{1❘{\sum\limits_{s \in {\{{{- 1},1}\}}}{\int{{\psi_{i,s}^{*}(r)}{\psi_{j,s}(r)}dr}}}} = \delta_{ij}}} \right\},} & (42) \end{matrix}$

where for each function ψ_(i,s) the index i labels the spatial degree of freedom and s labels the spin degree of freedom. Here the integer N is assumed to be even. The set of functions in Eqn. 42 is called a basis set and each element of the set is a spin orbital. For simplifying the notation from here on, we will use i, j, etc., to index combined spatial and spin degrees of freedom. The basis set in Eqn. 102 then becomes {ψ_(i)}. for i=1 to N. The introduction of a basis set allows one to discretize the problem. By introducing ladder operators {a_(i), a_(i) ^(†), i=1, . . . , N} satisfying fermionic commutation relationship {a_(i), a_(j) ^(†)}≡a_(i)a_(j) ^(†)+a_(j) ^(†)a′_(i)=δ_(ij)I, {a_(i) ^(†), a_(j) ^(†)}={a_(i), a_(j)}=0, {a_(i) ^(†), a_(j) ^(†)}{a_(i), a_(j)}=0, we can transform the Hamiltonian in Eqn. 1 into a second-quantized form:

$\begin{matrix} {{H = {{\sum\limits_{pq}{h_{pq}a_{p}^{\dagger}a_{q}}} + {\sum\limits_{pqrs}{h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}}}}},} & (43) \end{matrix}$

where the coefficients h_(pq) and h_(pq), are one-electron and two-electron integrals which can be computed efficiently classically.

In one construction, the classical computer 402 of FIG. 4 includes a first Hamiltonian 406, also referred to herein as an “electronic structure Hamiltonian.” The memory of the classical computer 402 also include data 408 representing a basis set, defined by Eqn. 42 above. The classical computer 402 transforms 410 the first Hamiltonian 406 into a second-quantized Hamiltonian 412, as shown in Eqn. 43 above.

The classical computer 402 then transforms the Hamiltonian 412 into a form that is more amenable for measurements on the quantum computer 404. The classical computer 402 does this by transforming 414 the Hamiltonian 412 into a linear combination of tensor products of Pauli operators:

$\begin{matrix} {{H = {{{\sum\limits_{pq}{h_{pq}O_{pq}}} + {\sum\limits_{pqrs}{h_{pqrs}O_{pqrs}}}} = {\sum\limits_{i}{h_{i}H_{i}}}}},} & (44) \end{matrix}$

where each operator O_(pq) and O_(pqrs) are sums of N-fold tensor products of operators from {I,X,Y,Z} that are mapped from a_(p) ^(†)a_(q) and a_(p) ^(†)a_(q) ^(†)a_(r)a_(s), respectively, under the isospectral transformation. In the second equality, we use an index i to go over all such n-fold tensor products of Pauli operators H_(i), representing the Hamiltonian H as a weighted sum of such terms. The result is a qubit-mapped Hamiltonian 416.

The classical computer 402 generates 418 a reference state 420 based on the first Hamiltonian 406 and the basis set 408. The classical computer 402 prepares 422 a reference state 424 on the quantum computer 404 based on the reference state 420 on the classical computer 402. Optionally, output 440 from measured energy 440 may also be provided as an input to reference state 420, as illustrated in FIG. 4.

The classical computer 402 generates coupled cluster amplitudes 426 (such as by using many-body perturbation theory, as shown in FIG. 4), which form the set of parameters that are tuned by the classical algorithms executed by the classical computer to minimize the energy expectation with respect to the Hamiltonian 416. The classical computer 402 generates 428 a quantum circuit 430 on the quantum computer 404 using the coupled cluster amplitudes 426 as parameters, such as by using Eqn. 48 (see below). The form of the operators O_(i) in Eqn. 47 depends on the isospectral mapping 432. In another construction, the isospectral mapping 432 is disposed in the location of arrow 414 between second-quantized Hamiltonian 412 and qubit-mapped Hamiltonian 416.

The quantum computer 404 executes 434 the quantum circuit 430 to prepare an ansatz state 436, such as by using Eqn. 49 (see below). The quantum computer 404 measures 438 the energy of the ansatz state 436, thereby producing a quantum energy measurement 440, such as by using Eqn. 45. The classical computer 402 converts 442 the quantum energy measurement 440 into a form suitable for storage on the classical computer 402, where the energy measurement takes the form of output energy 444 upon convergence of the optimization algorithm 442.

The above construction depicted in FIG. 4 may also be expressed as a framework of variational quantum eigensolver (VQE) using noisy intermediate scale quantum devices in conjunction with a classical computer for approximately preparing the ground state of quantum Hamiltonians (of the form in Eqn. 44 above) as well as approximately measuring the ground state energy.

The basic idea is to use a parametrized quantum circuit U(θ) to prepare an ansatz state |ψ_(ansatz)

=U(θ)|Φ_(ref)

from a reference state |Φ_(ref)

. To approximate the ground state of a Hamiltonian H, one uses a classical computer to optimize θ such that it minimizes energy of H:

ψ_(ansatz)(θ)|H|ψ _(ansatz)(θ)

=Σ_(i) h _(i)

ψ_(ansatz)(θ)|H _(i)|ψ_(ansatz)(θ)

  (45)

measured with respect to the ansatz state.

A reference state |Φ_(ref)

is suitable to be prepared on a quantum computer because it is a computational basis state (i.e., a product state where each qubit is either in |0

or |1

). The choice for ansatz state should be one which is likely beyond efficient classical computation, such as unitary coupled cluster: |ψ_(uCC)

=exp(T−T^(†))|Φ_(ref)

. Here, the cluster operator T is defined with a preexisting partition of spin orbitals in the reference state into those that are occupied by an electron (or equivalently the corresponding qubit is |1

, the set of which is labeled as occ) and those that are not occupied by an electron (or equivalently the corresponding qubit is |0

, the set of which is labeled as virt for “virtual orbitals”). The operator T is then written as a weighted combination of different levels of excitations from occupied spin orbitals to virtual spin orbitals:

$\begin{matrix} {T = {{T_{1} + T_{2} + \ldots} = {{\sum\limits_{{i \in {occ}},{a \in {virt}}}{t_{i}^{a}a_{a}^{\dagger}a_{i}}} + {\sum\limits_{i,{j \in {occ}},a,{b \in {virt}}}{t_{ij}^{ab}a_{a}^{\dagger}a_{b}^{\dagger}a_{i}a_{j}}} + {\ldots\mspace{14mu}.}}}} & (46) \end{matrix}$

Here, T₁ captures all single excitations (i.e., excitations of a single electron from spin orbital i to a), T₂ captures all double excitations (i.e., excitations of two electrons from spin orbitals i and j to a and b), etc. The coefficients t_(i) ^(a) and t_(ij) ^(ab) are coupled cluster amplitudes for single and double excitations, respectively. In both the uCC setting and the setting using vCC ansatz that is discussed here, these amplitudes form the set θ of parameters that are tuned by the classical algorithms to minimize the energy expectation with respect to the Hamiltonian.

The isospectral transformation which maps the Hamiltonian H from Eqns. 43 and 44 above can also be applied to the cluster operator T in Eqn. 46, producing an operator of the form:

$\begin{matrix} {T = {{{\sum\limits_{{i \in {occ}},{a \in {virt}}}{t_{i}^{a}O_{i}^{a}}} + {\sum\limits_{i,{j \in {occ}},a,{b \in {{vi}rt}}}{r_{ij}^{ab}O_{ij}^{ab}}} + \ldots} = {\sum\limits_{i}{t_{i}O_{i}}}}} & (47) \end{matrix}$

where each operator O_(ai) and O_(abij) has the same definition as in Eqn. 44. Because the set occ of occupied orbitals and the set virt of virtual orbitals are by definition disjoint and the isospectral transformation preserves the fermionic commutation relationship, all of the operators in {O_(i) ^(a),O_(ij) ^(ab), . . . } pairwise commute with each other. In the second equality in Eqn. 47 groups all of the indexing in {O_(i) ^(a),O_(ij) ^(ab), . . . } under one index i for simplifying notation. A crucial distinction between the grouping here from Eqn. 44 is that while each H_(i) in Eqn. 44 is an N-fold tensor product of operators from {I,X,Y,Z}, each O_(i) is an element of the set {O_(i) ^(a),O_(ij) ^(ab), . . . }, which could be a linear combination of N-fold tensor products of operators from {I,X,Y,Z}.

Since all of the O_(i) operators commute pairwise, and O_(i) ²=0, it is observed that for any positive integer r

$\begin{matrix} {{\exp(T)} = {\left\lbrack {\exp\left( {\frac{1}{r}T} \right)} \right\rbrack^{r} = {\left\lbrack {\prod\limits_{i}{\exp\left( {\frac{t_{i}}{r}O_{i}} \right)}} \right\rbrack^{r} = {\left\lbrack {\prod\limits_{i}\left( {I + {\frac{t_{i}}{r}O_{i}}} \right)} \right\rbrack^{r}.}}}} & (108) \end{matrix}$

Hence, to prepare the vCC ansatz state

$\begin{matrix} {\left. {❘\psi_{vCC}} \right\rangle = \frac{\left. {{\exp(T)}❘\Phi_{ref}} \right\rangle}{\left. {{\exp(T)}❘\Phi_{ref}} \right\rangle }} & (109) \end{matrix}$

one needs to consider implementation of each non-unitary operation

$I + {\frac{t_{i}}{r}O_{i}}$

in Eqn. 48. The parameter r can be chosen such that

$I + {\frac{t_{i}}{r}O_{i}}$

is arbitrarily close to identity. Then

$I + {\frac{t_{i}}{r}O_{i}}$

could be realized probabilistically using well established techniques for approximately implementing linear combination of unitaries (LCU).

Referring to FIG. 5, a flowchart is shown of a method 500 performed by the system 400 of FIG. 4 according to one embodiment of the present invention. In general, the method implements, on the quantum computer 404, a non-unitary operation of the form I+αU, where I is the identity operator, a is a scalar, and U is a unitary operator. The quantum computer 404 includes a plurality of qubits, including an ancilla qubit. The classical computer 402 includes a processor, a non-transitory computer-readable medium, and computer program instructions stored in the non-transitory computer-readable medium. The computer program instructions are executable by the processor to perform the method 500. The method includes: generating and storing, in the non-transitory computer-readable medium of the classical computer 402, computer-readable data representing a description of a first quantum circuit W which, when executed by the quantum computer 404, probabilistically realizes the non-unitary operation by the technique of linear combination of unitaries (FIG. 5, operation 502); generating and storing, in the non-transitory computer-readable medium of the classical computer 402, computer-readable data representing a description of a second quantum circuit, the second quantum circuit comprising a sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1) parametrized by an integer d, wherein S₀=WRW^(†)RW and wherein

$R = \begin{pmatrix} 1 & 0 \\ 0 & e^{i{\pi/3}} \end{pmatrix}$

is a single-qubit rotation applied to the ancilla qubit (FIG. 5, operation 504).

The method 500 may further include: on the quantum computer 404, executing the first quantum circuit to probabilistically realize the non-unitary operation by the technique of linear combination of unitaries (FIG. 5, operation 506). The method may further include: on the quantum computer 404, executing the second quantum circuit, comprising executing the sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1) (FIG. 5, operation 508).

It is to be understood that although the invention has been described above in terms of particular embodiments, the foregoing embodiments are provided as illustrative only, and do not limit or define the scope of the invention. Various other embodiments, including but not limited to the following, are also within the scope of the claims. For example, elements and components described herein may be further divided into additional components or joined together to form fewer components for performing the same functions.

Various physical embodiments of a quantum computer are suitable for use according to the present disclosure. In general, the fundamental data storage unit in quantum computing is the quantum bit, or qubit. The qubit is a quantum-computing analog of a classical digital computer system bit. A classical bit is considered to occupy, at any given point in time, one of two possible states corresponding to the binary digits (bits) 0 or 1. By contrast, a qubit is implemented in hardware by a physical medium with quantum-mechanical characteristics. Such a medium, which physically instantiates a qubit, may be referred to herein as a “physical instantiation of a qubit,” a “physical embodiment of a qubit,” a “medium embodying a qubit,” or similar terms, or simply as a “qubit,” for ease of explanation. It should be understood, therefore, that references herein to “qubits” within descriptions of embodiments of the present invention refer to physical media which embody qubits.

Each qubit has an infinite number of different potential quantum-mechanical states. When the state of a qubit is physically measured, the measurement produces one of two different basis states resolved from the state of the qubit. Thus, a single qubit can represent a one, a zero, or any quantum superposition of those two qubit states; a pair of qubits can be in any quantum superposition of 4 orthogonal basis states; and three qubits can be in any superposition of 8 orthogonal basis states. The function that defines the quantum-mechanical states of a qubit is known as its wavefunction. The wavefunction also specifies the probability distribution of outcomes for a given measurement. A qubit, which has a quantum state of dimension two (i.e., has two orthogonal basis states), may be generalized to a d-dimensional “qudit,” where d may be any integral value, such as 2, 3, 4, or higher. In the general case of a qudit, measurement of the qudit produces one of d different basis states resolved from the state of the qudit. Any reference herein to a qubit should be understood to refer more generally to an d-dimensional qudit with any value of d.

Although certain descriptions of qubits herein may describe such qubits in terms of their mathematical properties, each such qubit may be implemented in a physical medium in any of a variety of different ways. Examples of such physical media include superconducting material, trapped ions, photons, optical cavities, individual electrons trapped within quantum dots, point defects in solids (e.g., phosphorus donors in silicon or nitrogen-vacancy centers in diamond), molecules (e.g., alanine, vanadium complexes), or aggregations of any of the foregoing that exhibit qubit behavior, that is, comprising quantum states and transitions therebetween that can be controllably induced or detected.

For any given medium that implements a qubit, any of a variety of properties of that medium may be chosen to implement the qubit. For example, if electrons are chosen to implement qubits, then the x component of its spin degree of freedom may be chosen as the property of such electrons to represent the states of such qubits. Alternatively, the y component, or the z component of the spin degree of freedom may be chosen as the property of such electrons to represent the state of such qubits. This is merely a specific example of the general feature that for any physical medium that is chosen to implement qubits, there may be multiple physical degrees of freedom (e.g., the x, y, and z components in the electron spin example) that may be chosen to represent 0 and 1. For any particular degree of freedom, the physical medium may controllably be put in a state of superposition, and measurements may then be taken in the chosen degree of freedom to obtain readouts of qubit values.

Certain implementations of quantum computers, referred as gate model quantum computers, comprise quantum gates. In contrast to classical gates, there is an infinite number of possible single-qubit quantum gates that change the state vector of a qubit. Changing the state of a qubit state vector typically is referred to as a single-qubit rotation, and may also be referred to herein as a state change or a single-qubit quantum-gate operation. A rotation, state change, or single-qubit quantum-gate operation may be represented mathematically by a unitary 2×2 matrix with complex elements. A rotation corresponds to a rotation of a qubit state within its Hilbert space, which may be conceptualized as a rotation of the Bloch sphere. (As is well-known to those having ordinary skill in the art, the Bloch sphere is a geometrical representation of the space of pure states of a qubit.) Multi-qubit gates alter the quantum state of a set of qubits. For example, two-qubit gates rotate the state of two qubits as a rotation in the four-dimensional Hilbert space of the two qubits. (As is well-known to those having ordinary skill in the art, a Hilbert space is an abstract vector space possessing the structure of an inner product that allows length and angle to be measured. Furthermore, Hilbert spaces are complete: there are enough limits in the space to allow the techniques of calculus to be used.)

A quantum circuit may be specified as a sequence of quantum gates. As described in more detail below, the term “quantum gate,” as used herein, refers to the application of a gate control signal (defined below) to one or more qubits to cause those qubits to undergo certain physical transformations and thereby to implement a logical gate operation. To conceptualize a quantum circuit, the matrices corresponding to the component quantum gates may be multiplied together in the order specified by the gate sequence to produce a 2n×2n complex matrix representing the same overall state change on n qubits. A quantum circuit may thus be expressed as a single resultant operator. However, designing a quantum circuit in terms of constituent gates allows the design to conform to a standard set of gates, and thus enable greater ease of deployment. A quantum circuit thus corresponds to a design for actions taken upon the physical components of a quantum computer.

A given variational quantum circuit may be parameterized in a suitable device-specific manner. More generally, the quantum gates making up a quantum circuit may have an associated plurality of tuning parameters. For example, in embodiments based on optical switching, tuning parameters may correspond to the angles of individual optical elements.

In certain embodiments of quantum circuits, the quantum circuit includes both one or more gates and one or more measurement operations. Quantum computers implemented using such quantum circuits are referred to herein as implementing “measurement feedback.” For example, a quantum computer implementing measurement feedback may execute the gates in a quantum circuit and then measure only a subset (i.e., fewer than all) of the qubits in the quantum computer, and then decide which gate(s) to execute next based on the outcome(s) of the measurement(s). In particular, the measurement(s) may indicate a degree of error in the gate operation(s), and the quantum computer may decide which gate(s) to execute next based on the degree of error. The quantum computer may then execute the gate(s) indicated by the decision. This process of executing gates, measuring a subset of the qubits, and then deciding which gate(s) to execute next may be repeated any number of times. Measurement feedback may be useful for performing quantum error correction, but is not limited to use in performing quantum error correction. For every quantum circuit, there is an error-corrected implementation of the circuit with or without measurement feedback.

Some embodiments described herein generate, measure, or utilize quantum states that approximate a target quantum state (e.g., a ground state of a Hamiltonian). As will be appreciated by those trained in the art, there are many ways to quantify how well a first quantum state “approximates” a second quantum state. In the following description, any concept or definition of approximation known in the art may be used without departing from the scope hereof. For example, when the first and second quantum states are represented as first and second vectors, respectively, the first quantum state approximates the second quantum state when an inner product between the first and second vectors (called the “fidelity” between the two quantum states) is greater than a predefined amount (typically labeled E). In this example, the fidelity quantifies how “close” or “similar” the first and second quantum states are to each other. The fidelity represents a probability that a measurement of the first quantum state will give the same result as if the measurement were performed on the second quantum state. Proximity between quantum states can also be quantified with a distance measure, such as a Euclidean norm, a Hamming distance, or another type of norm known in the art. Proximity between quantum states can also be defined in computational terms. For example, the first quantum state approximates the second quantum state when a polynomial time-sampling of the first quantum state gives some desired information or property that it shares with the second quantum state.

Not all quantum computers are gate model quantum computers. Embodiments of the present invention are not limited to being implemented using gate model quantum computers. As an alternative example, embodiments of the present invention may be implemented, in whole or in part, using a quantum computer that is implemented using a quantum annealing architecture, which is an alternative to the gate model quantum computing architecture. More specifically, quantum annealing (QA) is a metaheuristic for finding the global minimum of a given objective function over a given set of candidate solutions (candidate states), by a process using quantum fluctuations.

FIG. X2 shows a diagram illustrating operations typically performed by a computer system X50 which implements quantum annealing. The system X50 includes both a quantum computer X52 and a classical computer X54. Operations shown on the left of the dashed vertical line X56 typically are performed by the quantum computer X52, while operations shown on the right of the dashed vertical line X56 typically are performed by the classical computer X54.

Quantum annealing starts with the classical computer X54 generating an initial Hamiltonian X60 and a final Hamiltonian X62 based on a computational problem X58 to be solved, and providing the initial Hamiltonian X60, the final Hamiltonian X62 and an annealing schedule X70 as input to the quantum computer X52. The quantum computer X52 prepares a well-known initial state X66 (FIG. X2, operation X64), such as a quantum-mechanical superposition of all possible states (candidate states) with equal weights, based on the initial Hamiltonian X60. The classical computer X54 provides the initial Hamiltonian X60, a final Hamiltonian X62, and an annealing schedule X70 to the quantum computer X52. The quantum computer X52 starts in the initial state X66, and evolves its state according to the annealing schedule X70 following the time-dependent Schrödinger equation, a natural quantum-mechanical evolution of physical systems (FIG. X2, operation X68). More specifically, the state of the quantum computer X52 undergoes time evolution under a time-dependent Hamiltonian, which starts from the initial Hamiltonian X60 and terminates at the final Hamiltonian X62. If the rate of change of the system Hamiltonian is slow enough, the system stays close to the ground state of the instantaneous Hamiltonian. If the rate of change of the system Hamiltonian is accelerated, the system may leave the ground state temporarily but produce a higher likelihood of concluding in the ground state of the final problem Hamiltonian, i.e., diabatic quantum computation. At the end of the time evolution, the set of qubits on the quantum annealer is in a final state X72, which is expected to be close to the ground state of the classical Ising model that corresponds to the solution to the original optimization problem X58. An experimental demonstration of the success of quantum annealing for random magnets was reported immediately after the initial theoretical proposal.

The final state X72 of the quantum computer X54 is measured, thereby producing results X76 (i.e., measurements) (FIG. X2, operation X74). The measurement operation X74 may be performed, for example, in any of the ways disclosed herein, such as in any of the ways disclosed herein in connection with the measurement unit W10 in FIG. W. The classical computer X54 performs postprocessing on the measurement results X76 to produce output X80 representing a solution to the original computational problem X58 (FIG. X2, operation X78).

As yet another alternative example, embodiments of the present invention may be implemented, in whole or in part, using a quantum computer that is implemented using a one-way quantum computing architecture, also referred to as a measurement-based quantum computing architecture, which is another alternative to the gate model quantum computing architecture. More specifically, the one-way or measurement based quantum computer (MBQC) is a method of quantum computing that first prepares an entangled resource state, usually a cluster state or graph state, then performs single qubit measurements on it. It is “one-way” because the resource state is destroyed by the measurements.

The outcome of each individual measurement is random, but they are related in such a way that the computation always succeeds. In general the choices of basis for later measurements need to depend on the results of earlier measurements, and hence the measurements cannot all be performed at the same time.

Any of the functions disclosed herein may be implemented using means for performing those functions. Such means include, but are not limited to, any of the components disclosed herein, such as the computer-related components described below.

Referring to FIG. W, a diagram is shown of a system W00 implemented according to one embodiment of the present invention. Referring to FIG. X1, a flowchart is shown of a method X00 performed by the system W00 of FIG. W according to one embodiment of the present invention. The system W00 includes a quantum computer W02. The quantum computer W02 includes a plurality of qubits W04, which may be implemented in any of the ways disclosed herein. There may be any number of qubits W04 in the quantum computer W04. For example, the qubits W04 may include or consist of no more than 2 qubits, no more than 4 qubits, no more than 8 qubits, no more than 16 qubits, no more than 32 qubits, no more than 64 qubits, no more than 128 qubits, no more than 256 qubits, no more than 512 qubits, no more than 1024 qubits, no more than 2048 qubits, no more than 4096 qubits, or no more than 8192 qubits. These are merely examples, in practice there may be any number of qubits W04 in the quantum computer W02.

There may be any number of gates in a quantum circuit. However, in some embodiments the number of gates may be at least proportional to the number of qubits W04 in the quantum computer W02. In some embodiments the gate depth may be no greater than the number of qubits W04 in the quantum computer W02, or no greater than some linear multiple of the number of qubits W04 in the quantum computer W02 (e.g., 2, 3, 4, 5, 6, or 7).

The qubits W04 may be interconnected in any graph pattern. For example, they be connected in a linear chain, a two-dimensional grid, an all-to-all connection, any combination thereof, or any subgraph of any of the preceding.

As will become clear from the description below, although element W02 is referred to herein as a “quantum computer,” this does not imply that all components of the quantum computer W02 leverage quantum phenomena. One or more components of the quantum computer W02 may, for example, be classical (i.e., non-quantum components) components which do not leverage quantum phenomena.

The quantum computer W02 includes a control unit W06, which may include any of a variety of circuitry and/or other machinery for performing the functions disclosed herein. The control unit W06 may, for example, consist entirely of classical components. The control unit W06 generates and provides as output one or more control signals W08 to the qubits W04. The control signals W08 may take any of a variety of forms, such as any kind of electromagnetic signals, such as electrical signals, magnetic signals, optical signals (e.g., laser pulses), or any combination thereof.

For example:

In embodiments in which some or all of the qubits W04 are implemented as photons (also referred to as a “quantum optical” implementation) that travel along waveguides, the control unit W06 may be a beam splitter (e.g., a heater or a mirror), the control signals W08 may be signals that control the heater or the rotation of the mirror, the measurement unit W10 may be a photodetector, and the measurement signals W12 may be photons.

In embodiments in which some or all of the qubits W04 are implemented as charge type qubits (e.g., transmon, X-mon, G-mon) or flux-type qubits (e.g., flux qubits, capacitively shunted flux qubits) (also referred to as a “circuit quantum electrodynamic” (circuit QED) implementation), the control unit W06 may be a bus resonator activated by a drive, the control signals W08 may be cavity modes, the measurement unit W10 may be a second resonator (e.g., a low-Q resonator), and the measurement signals W12 may be voltages measured from the second resonator using dispersive readout techniques.

In embodiments in which some or all of the qubits W04 are implemented as superconducting circuits, the control unit W06 may be a circuit QED-assisted control unit or a direct capacitive coupling control unit or an inductive capacitive coupling control unit, the control signals W08 may be cavity modes, the measurement unit W10 may be a second resonator (e.g., a low-Q resonator), and the measurement signals W12 may be voltages measured from the second resonator using dispersive readout techniques.

In embodiments in which some or all of the qubits W04 are implemented as trapped ions (e.g., electronic states of, e.g., magnesium ions), the control unit W06 may be a laser, the control signals W08 may be laser pulses, the measurement unit W10 may be a laser and either a CCD or a photodetector (e.g., a photomultiplier tube), and the measurement signals W12 may be photons.

In embodiments in which some or all of the qubits W04 are implemented using nuclear magnetic resonance (NMR) (in which case the qubits may be molecules, e.g., in liquid or solid form), the control unit W06 may be a radio frequency (RF) antenna, the control signals W08 may be RF fields emitted by the RF antenna, the measurement unit W10 may be another RF antenna, and the measurement signals W12 may be RF fields measured by the second RF antenna.

In embodiments in which some or all of the qubits W04 are implemented as nitrogen-vacancy centers (NV centers), the control unit W06 may, for example, be a laser, a microwave antenna, or a coil, the control signals W08 may be visible light, a microwave signal, or a constant electromagnetic field, the measurement unit W10 may be a photodetector, and the measurement signals W12 may be photons.

In embodiments in which some or all of the qubits W04 are implemented as two-dimensional quasiparticles called “anyons” (also referred to as a “topological quantum computer” implementation), the control unit W06 may be nanowires, the control signals W08 may be local electrical fields or microwave pulses, the measurement unit W10 may be superconducting circuits, and the measurement signals W12 may be voltages.

In embodiments in which some or all of the qubits W04 are implemented as semiconducting material (e.g., nanowires), the control unit W06 may be microfabricated gates, the control signals W08 may be RF or microwave signals, the measurement unit W10 may be microfabricated gates, and the measurement signals W12 may be RF or microwave signals.

Although not shown explicitly in FIG. W and not required, the measurement unit W10 may provide one or more feedback signals W14 to the control unit W06 based on the measurement signals W12. For example, quantum computers referred to as “one-way quantum computers” or “measurement-based quantum computers” utilize such feedback W14 from the measurement unit W10 to the control unit W06. Such feedback W14 is also necessary for the operation of fault-tolerant quantum computing and error correction.

The control signals W08 may, for example, include one or more state preparation signals which, when received by the qubits W04, cause some or all of the qubits W04 to change their states. Such state preparation signals constitute a quantum circuit also referred to as an “ansatz circuit.” The resulting state of the qubits W04 is referred to herein as an “initial state” or an “ansatz state.” The process of outputting the state preparation signal(s) to cause the qubits W04 to be in their initial state is referred to herein as “state preparation” (FIG. X1, section X06). A special case of state preparation is “initialization,” also referred to as a “reset operation,” in which the initial state is one in which some or all of the qubits W04 are in the “zero” state i.e. the default single-qubit state. More generally, state preparation may involve using the state preparation signals to cause some or all of the qubits W04 to be in any distribution of desired states. In some embodiments, the control unit W06 may first perform initialization on the qubits W04 and then perform preparation on the qubits W04, by first outputting a first set of state preparation signals to initialize the qubits W04, and by then outputting a second set of state preparation signals to put the qubits W04 partially or entirely into non-zero states.

Another example of control signals W08 that may be output by the control unit W06 and received by the qubits W04 are gate control signals. The control unit W06 may output such gate control signals, thereby applying one or more gates to the qubits W04. Applying a gate to one or more qubits causes the set of qubits to undergo a physical state change which embodies a corresponding logical gate operation (e.g., single-qubit rotation, two-qubit entangling gate or multi-qubit operation) specified by the received gate control signal. As this implies, in response to receiving the gate control signals, the qubits W04 undergo physical transformations which cause the qubits W04 to change state in such a way that the states of the qubits W04, when measured (see below), represent the results of performing logical gate operations specified by the gate control signals. The term “quantum gate,” as used herein, refers to the application of a gate control signal to one or more qubits to cause those qubits to undergo the physical transformations described above and thereby to implement a logical gate operation.

It should be understood that the dividing line between state preparation (and the corresponding state preparation signals) and the application of gates (and the corresponding gate control signals) may be chosen arbitrarily. For example, some or all the components and operations that are illustrated in FIGS. W and X as elements of “state preparation” may instead be characterized as elements of gate application.

Conversely, for example, some or all of the components and operations that are illustrated in FIGS. W and X as elements of “gate application” may instead be characterized as elements of state preparation. As one particular example, the system and method of FIGS. W and X may be characterized as solely performing state preparation followed by measurement, without any gate application, where the elements that are described herein as being part of gate application are instead considered to be part of state preparation. Conversely, for example, the system and method of FIGS. W and X may be characterized as solely performing gate application followed by measurement, without any state preparation, and where the elements that are described herein as being part of state preparation are instead considered to be part of gate application.

The quantum computer W02 also includes a measurement unit W10, which performs one or more measurement operations on the qubits W04 to read out measurement signals W12 (also referred to herein as “measurement results”) from the qubits W04, where the measurement results W12 are signals representing the states of some or all of the qubits W04. In practice, the control unit

W06 and the measurement unit W10 may be entirely distinct from each other, or contain some components in common with each other, or be implemented using a single unit (i.e., a single unit may implement both the control unit W06 and the measurement unit W10). For example, a laser unit may be used both to generate the control signals W08 and to provide stimulus (e.g., one or more laser beams) to the qubits W04 to cause the measurement signals W12 to be generated.

In general, the quantum computer W02 may perform various operations described above any number of times. For example, the control unit W06 may generate one or more control signals W08, thereby causing the qubits W04 to perform one or more quantum gate operations. The measurement unit W10 may then perform one or more measurement operations on the qubits W04 to read out a set of one or more measurement signals W12. The measurement unit W10 may repeat such measurement operations on the qubits W04 before the control unit W06 generates additional control signals W08, thereby causing the measurement unit W10 to read out additional measurement signals W12 resulting from the same gate operations that were performed before reading out the previous measurement signals W12. The measurement unit W10 may repeat this process any number of times to generate any number of measurement signals W12 corresponding to the same gate operations. The quantum computer W02 may then aggregate such multiple measurements of the same gate operations in any of a variety of ways.

After the measurement unit W10 has performed one or more measurement operations on the qubits W04 after they have performed one set of gate operations, the control unit W06 may generate one or more additional control signals W08, which may differ from the previous control signals W08, thereby causing the qubits W04 to perform one or more additional quantum gate operations, which may differ from the previous set of quantum gate operations. The process described above may then be repeated, with the measurement unit W10 performing one or more measurement operations on the qubits W04 in their new states (resulting from the most recently-performed gate operations).

In general, the system W00 may implement a plurality of quantum circuits as follows. For each quantum circuit C in the plurality of quantum circuits (FIG. X1, operation X02), the system W00 performs a plurality of “shots” on the qubits W04. The meaning of a shot will become clear from the description that follows. For each shot S in the plurality of shots (FIG. X1, operation X04), the system W00 prepares the state of the qubits W04 (FIG. X1, section X06). More specifically, for each quantum gate G in quantum circuit C (FIG. X1, operation X10), the system W00 applies quantum gate G to the qubits W04 (FIG. X1, operations X12 and X14).

Then, for each of the qubits Q W04 (FIG. X1, operation X16), the system W00 measures the qubit Q to produce measurement output representing a current state of qubit Q (FIG. X1, operations X18 and X20).

The operations described above are repeated for each shot S (FIG. X1, operation X22), and circuit C (FIG. X1, operation X24). As the description above implies, a single “shot” involves preparing the state of the qubits W04 and applying all of the quantum gates in a circuit to the qubits W04 and then measuring the states of the qubits W04; and the system W00 may perform multiple shots for one or more circuits.

Referring to FIG. Y, a diagram is shown of a hybrid classical quantum computer (HQC) Y00 implemented according to one embodiment of the present invention. The HQC Y00 includes a quantum computer component W02 (which may, for example, be implemented in the manner shown and described in connection with FIG. W) and a classical computer component Y06. The classical computer component may be a machine implemented according to the general computing model established by John Von Neumann, in which programs are written in the form of ordered lists of instructions and stored within a classical (e.g., digital) memory Y10 and executed by a classical (e.g., digital) processor Y08 of the classical computer. The memory Y10 is classical in the sense that it stores data in a storage medium in the form of bits, which have a single definite binary state at any point in time. The bits stored in the memory Y10 may, for example, represent a computer program. The classical computer component Y04 typically includes a bus Y14. The processor Y08 may read bits from and write bits to the memory Y10 over the bus Y14. For example, the processor Y08 may read instructions from the computer program in the memory Y10, and may optionally receive input data Y16 from a source external to the computer Y02, such as from a user input device such as a mouse, keyboard, or any other input device. The processor Y08 may use instructions that have been read from the memory Y10 to perform computations on data read from the memory Y10 and/or the input Y16, and generate output from those instructions. The processor Y08 may store that output back into the memory Y10 and/or provide the output externally as output data Y18 via an output device, such as a monitor, speaker, or network device.

The quantum computer component W02 may include a plurality of qubits W04, as described above in connection with FIG. W. A single qubit may represent a one, a zero, or any quantum superposition of those two qubit states. The classical computer component Y04 may provide classical state preparation signals Y32 to the quantum computer W02, in response to which the quantum computer W02 may prepare the states of the qubits W04 in any of the ways disclosed herein, such as in any of the ways disclosed in connection with FIGS. W and X.

Once the qubits W04 have been prepared, the classical processor Y08 may provide classical control signals Y34 to the quantum computer W02, in response to which the quantum computer W02 may apply the gate operations specified by the control signals Y32 to the qubits W04, as a result of which the qubits W04 arrive at a final state. The measurement unit W10 in the quantum computer W02 (which may be implemented as described above in connection with FIGS. W and X) may measure the states of the qubits W04 and produce measurement output Y38 representing the collapse of the states of the qubits W04 into one of their eigenstates. As a result, the measurement output Y38 includes or consists of bits and therefore represents a classical state. The quantum computer W02 provides the measurement output Y38 to the classical processor Y08. The classical processor Y08 may store data representing the measurement output Y38 and/or data derived therefrom in the classical memory Y10.

The steps described above may be repeated any number of times, with what is described above as the final state of the qubits W04 serving as the initial state of the next iteration. In this way, the classical computer Y04 and the quantum computer W02 may cooperate as co-processors to perform joint computations as a single computer system.

Although certain functions may be described herein as being performed by a classical computer and other functions may be described herein as being performed by a quantum computer, these are merely examples and do not constitute limitations of the present invention. A subset of the functions which are disclosed herein as being performed by a quantum computer may instead be performed by a classical computer. For example, a classical computer may execute functionality for emulating a quantum computer and provide a subset of the functionality described herein, albeit with functionality limited by the exponential scaling of the simulation. Functions which are disclosed herein as being performed by a classical computer may instead be performed by a quantum computer.

The techniques described above may be implemented, for example, in hardware, in one or more computer programs tangibly stored on one or more computer-readable media, firmware, or any combination thereof, such as solely on a quantum computer, solely on a classical computer, or on a hybrid classical quantum (HQC) computer. The techniques disclosed herein may, for example, be implemented solely on a classical computer, in which the classical computer emulates the quantum computer functions disclosed herein.

The techniques described above may be implemented in one or more computer programs executing on (or executable by) a programmable computer (such as a classical computer, a quantum computer, or an HQC) including any combination of any number of the following: a processor, a storage medium readable and/or writable by the processor (including, for example, volatile and non-volatile memory and/or storage elements), an input device, and an output device. Program code may be applied to input entered using the input device to perform the functions described and to generate output using the output device.

Embodiments of the present invention include features which are only possible and/or feasible to implement with the use of one or more computers, computer processors, and/or other elements of a computer system. Such features are either impossible or impractical to implement mentally and/or manually. For example, the quantum circuit execution described herein as being performed by a quantum computer, which performs quantum operations on quantum states, which for moderately large systems (e.g., at least 50 qubits) would be infeasible or impossible to perform manually or even using a classical computer.

Any claims herein which affirmatively require a computer, a processor, a memory, or similar computer-related elements, are intended to require such elements, and should not be interpreted as if such elements are not present in or required by such claims. Such claims are not intended, and should not be interpreted, to cover methods and/or systems which lack the recited computer-related elements. For example, any method claim herein which recites that the claimed method is performed by a computer, a processor, a memory, and/or similar computer-related element, is intended to, and should only be interpreted to, encompass methods which are performed by the recited computer-related element(s). Such a method claim should not be interpreted, for example, to encompass a method that is performed mentally or by hand (e.g., using pencil and paper). Similarly, any product claim herein which recites that the claimed product includes a computer, a processor, a memory, and/or similar computer-related element, is intended to, and should only be interpreted to, encompass products which include the recited computer-related element(s). Such a product claim should not be interpreted, for example, to encompass a product that does not include the recited computer-related element(s).

In embodiments in which a classical computing component executes a computer program providing any subset of the functionality within the scope of the claims below, the computer program may be implemented in any programming language, such as assembly language, machine language, a high-level procedural programming language, or an object-oriented programming language. The programming language may, for example, be a compiled or interpreted programming language.

Each such computer program may be implemented in a computer program product tangibly embodied in a machine-readable storage device for execution by a computer processor, which may be either a classical processor or a quantum processor. Method steps of the invention may be performed by one or more computer processors executing a program tangibly embodied on a computer-readable medium to perform functions of the invention by operating on input and generating output. Suitable processors include, by way of example, both general and special purpose microprocessors. Generally, the processor receives (reads) instructions and data from a memory (such as a read-only memory and/or a random access memory) and writes (stores) instructions and data to the memory. Storage devices suitable for tangibly embodying computer program instructions and data include, for example, all forms of non-volatile memory, such as semiconductor memory devices, including EPROM, EEPROM, and flash memory devices; magnetic disks such as internal hard disks and removable disks; magneto-optical disks; and CD-ROMs. Any of the foregoing may be supplemented by, or incorporated in, specially-designed ASICs (application-specific integrated circuits) or FPGAs (Field-Programmable Gate Arrays). A classical computer can generally also receive (read) programs and data from, and write (store) programs and data to, a non-transitory computer-readable storage medium such as an internal disk (not shown) or a removable disk. These elements will also be found in a conventional desktop or workstation computer as well as other computers suitable for executing computer programs implementing the methods described herein, which may be used in conjunction with any digital print engine or marking engine, display monitor, or other raster output device capable of producing color or gray scale pixels on paper, film, display screen, or other output medium.

Any data disclosed herein may be implemented, for example, in one or more data structures tangibly stored on a non-transitory computer-readable medium (such as a classical computer-readable medium, a quantum computer-readable medium, or an HQC computer-readable medium). Embodiments of the invention may store such data in such data structure(s) and read such data from such data structure(s). 

1. A method performed by a classical computer for implementing, on a quantum computer, a non-unitary operation of the form I+αU, where I is the identity operator, α is a scalar, and U is a unitary operator, the quantum computer having a plurality of qubits, including an ancilla qubit, the classical computer including a processor, a non-transitory computer-readable medium, and computer program instructions stored in the non-transitory computer-readable medium, the computer program instructions being executable by the processor to perform the method, the method comprising: (A) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a first quantum circuit W which, when executed by the quantum computer, probabilistically realizes the non-unitary operation by the technique of linear combination of unitaries; (B) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a second quantum circuit, the second quantum circuit comprising a sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1) parametrized by an integer d, wherein S₀=WRW^(†)RW and wherein $R = \begin{pmatrix} 1 & 0 \\ 0 & e^{i{\pi/3}} \end{pmatrix}$  is a single-qubit rotation applied to the ancilla qubit.
 2. The method of claim 1, further comprising: (C) on the quantum computer, executing the first quantum circuit to probabilistically realize the non-unitary operation by the technique of linear combination of unitaries.
 3. The method of claim 2, further comprising: (D) on the quantum computer, executing the second quantum circuit, comprising executing the sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1).
 4. The method of claim 3, wherein (D) comprises, on the quantum computer: (D)(1) applying a single-qubit rotation R_(α) to transform the ancilla qubit into the state ${{\frac{1}{\sqrt{1 + \alpha}}\left. 0 \right\rangle} + {\sqrt{\frac{\alpha}{1 + \alpha}}\left. 1 \right\rangle}};$ (D)(2) applying controlled-U operator on a subset of the plurality of qubits, not including the ancilla qubit, conditioned on the ancilla qubit being in the state |1

; and (D)(3) applying R_(α) ^(†) on the ancilla qubit.
 5. The method of claim 1, further comprising generating and storing, in the non-transitory computer-readable medium, additional computer-readable data that, when executed on the quantum computer, causes the second quantum circuit to execute repeatedly, on the quantum computer, to perform a sequence of operations ${\prod_{i = 1}^{r}{\prod_{j = 1}^{k}\left( {1 + {\frac{t\kappa_{j}}{r}P_{j}}} \right)}},$ that approximates the operator e^(tT({right arrow over (κ)})) with t being a scalar and T({right arrow over (κ)})=Σ_(j=1) ^(k) κ_(j)P_(j) being the cluster operator which is a linear combination of operators P_(j) parametrized by weights {right arrow over (κ)}.
 6. The method of claim 5, further comprising: (C) performing mean-field approximation to generate and store, in the non-transitory computer-readable medium, computer-readable data representing a description of a third quantum circuit to prepare a reference state; (D) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a parametrized quantum circuit for approximating e^(tT({right arrow over (κ)})); (E) on the quantum computer, executing the third quantum circuit to prepare the reference state; (F) on the quantum computer, applying the parametrized quantum circuit to the reference state to generate the ansatz state |ψ({right arrow over (κ)})

; (G) using the classical computer in cooperation with the quantum computer to measure an energy of the ansatz |ψ({right arrow over (κ)})

; and (H) on the classical computer, iteratively tuning the parameters {right arrow over (κ)} is to minimize the energy of the ansatz |ψ({right arrow over (κ)})

.
 7. A system comprising: a classical computer the classical computer including a processor, a non-transitory computer-readable medium, and computer program instructions stored in the non-transitory computer-readable medium; a quantum computer comprising a plurality of qubits, including an ancilla qubit; wherein the computer program instructions, when executed by the processor, perform a method for implementing, on the quantum computer, a non-unitary operation of the form I+αU, where I is the identity operator, α is a scalar, and U is a unitary operator, the method comprising: (A) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a first quantum circuit W which, when executed by the quantum computer, probabilistically realizes the non-unitary operation by the technique of linear combination of unitaries; (B) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a second quantum circuit, the second quantum circuit comprising a sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1) parametrized by an integer d, wherein S₀=WRW^(†)RW and wherein $R = \begin{pmatrix} 1 & 0 \\ 0 & e^{i{\pi/3}} \end{pmatrix}$  is a single-qubit rotation applied to the ancilla qubit.
 8. The system of claim 7, wherein the method further comprises: (C) on the quantum computer, executing the first quantum circuit to probabilistically realize the non-unitary operation by the technique of linear combination of unitaries.
 9. The system of claim 8, wherein the method further comprises: (D) on the quantum computer, executing the second quantum circuit, comprising executing the sequence of quantum gates S_(d)=S_(d−1)RS_(d−1) ^(†)RS_(d−1).
 10. The system of claim 9, wherein (D) comprises, on the quantum computer: (D)(1) applying a single-qubit rotation R_(α) to transform the ancilla qubit into the state ${{\frac{1}{\sqrt{1 + \alpha}}\left. 0 \right\rangle} + {\sqrt{\frac{\alpha}{1 + \alpha}}\left. 1 \right\rangle}};$ (D)(2) applying controlled-U operator on a subset of the plurality of qubits, not including the ancilla qubit, conditioned on the ancilla qubit being in the state |1

; and (D)(3) applying R_(α) ^(†) on the ancilla qubit.
 11. The system of claim 7, wherein the method further comprises generating and storing, in the non-transitory computer-readable medium, additional computer-readable data that, when executed on the quantum computer, causes the second quantum circuit to execute repeatedly, on the quantum computer, to perform a sequence of operations ${\prod_{i = 1}^{r}{\prod_{j = 1}^{k}\left( {1 + {\frac{t\kappa_{j}}{r}P_{j}}} \right)}},$ that approximates the operator e^(tT({right arrow over (κ)})) with t being a scalar and T({right arrow over (κ)}))=Σ_(j=1) ^(k) κ_(j)P_(j) being the cluster operator which is a linear combination of operators P_(j) parametrized by weights {right arrow over (κ)}.
 12. The system of claim 11, wherein the method further comprises: (C) performing mean-field approximation to generate and store, in the non-transitory computer-readable medium, computer-readable data representing a description of a third quantum circuit to prepare a reference state; (D) generating and storing, in the non-transitory computer-readable medium, computer-readable data representing a description of a parametrized quantum circuit for approximating e^(tT({right arrow over (κ)})); (E) on the quantum computer, executing the third quantum circuit to prepare the reference state; (F) on the quantum computer, applying the parametrized quantum circuit to the reference state to generate the ansatz state |ψ({right arrow over (κ)})

; (G) using the classical computer in cooperation with the quantum computer to measure an energy of the ansatz |ψ({right arrow over (κ)})

; and (H) on the classical computer, iteratively tuning the parameters {right arrow over (K)} is to minimize the energy of the ansatz |ψ({right arrow over (κ)})

. 